#include "../include/SSS_coarsen.h"
//#include <__clang_cuda_device_functions.h>
#include <cstdlib>
#include <cuda_runtime.h>
#include <driver_types.h>
#include <numeric>
#include <iterator>
#include <iostream>
#include <vector>
#include <algorithm>
#include <bits/stdc++.h>
#include <metis.h>
#include <stdlib.h>






#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<stdint.h>
#include<stdarg.h>
#include<time.h>
#include<cuda_runtime.h>
#include<omp.h>
#include<sys/time.h>
#include<thrust/scan.h>
#include<thrust/sort.h>
#include<thrust/execution_policy.h>
#include<thrust/device_ptr.h>
#include<omp.h>
#include "bb_segsort.h"



extern double SoC_time ; 
extern double Coarsen_time ; 
extern double metis_time ;

#define block_size 64

int test_time;

clock_t t_contract_start,t_contract_end; 
double duration_contract=0;

clock_t t_match_start,t_match_end; 
double duration_match=0;

clock_t t_refine_start,t_refine_end; 
double duration_refine=0;

clock_t t_project_start,t_project_end; 
double duration_project=0;

clock_t t_coarsen_start,t_coarsen_end; 
double duration_coarsen=0;

clock_t t_uncoarsen_start,t_uncoarsen_end; 
double duration_uncoarsen=0;

clock_t t_init_start,t_init_end; 
double duration_init=0;

clock_t t_part_start,t_part_end; 
double duration_part=0;

//match
clock_t t_cmk1_start,t_cmk1_end; //cuda_match
double duration_cmk1=0;

clock_t t_cmk2_start,t_cmk2_end; //cuda_cleanv
double duration_cmk2=0;

clock_t t_cmk3_start,t_cmk3_end; //findc1
double duration_cmk3=0;

clock_t t_cmk4_start,t_cmk4_end; //findc2
double duration_cmk4=0;

clock_t t_cmk5_start,t_cmk5_end; //inclusive_scan
double duration_cmk5=0;

clock_t t_cmk6_start,t_cmk6_end; //findc2_5
double duration_cmk6=0;

clock_t t_cmk7_start,t_cmk7_end; //findc3
double duration_cmk7=0;

clock_t t_cmk8_start,t_cmk8_end; //findc4
double duration_cmk8=0;




clock_t t_cck1_start,t_cck1_end; //exclusive_scan
double duration_cck1=0;

clock_t t_cck2_start,t_cck2_end; //find_cnvtxsedge_original
double duration_cck2=0;

clock_t t_cck3_start,t_cck3_end; //bb_segsort
double duration_cck3=0;

clock_t t_cck4_start,t_cck4_end; //Sort_cnedges_part2
double duration_cck4=0;

clock_t t_cck5_start,t_cck5_end; //exclusive_scan
double duration_cck5=0;

clock_t t_cck6_start,t_cck6_end; //create_cadjncy
double duration_cck6=0;



clock_t t_rpk1_start,t_rpk1_end; //projectback
double duration_rpk1=0;


clock_t t_rrk1_start,t_rrk1_end; //find info
double duration_rrk1=0;

clock_t t_rrk2_start,t_rrk2_end; //find para
double duration_rrk2=0;

clock_t t_rrk3_start,t_rrk3_end; //exchange
double duration_rrk3=0;

clock_t t_gtoc_start,t_gtoc_end; //gputocpu
double duration_gtoc=0;

clock_t t_ctog_start,t_ctog_end; //cputogpu
double duration_ctog=0;


typedef struct ckrinfo_t {
 int id;              
 int ed;            	
 int nnbrs;          
 int inbr;            
} ckrinfo_t;

typedef struct vkrinfo_t {
 int nid;             
 int ned;            	
 int gv;            	
 int nnbrs;          	
 int inbr;            
} vkrinfo_t;

#define gk_max(a, b) ((a) >= (b) ? (a) : (b))
#define gk_min(a, b) ((a) >= (b) ? (b) : (a))
#define LTERM			(void **) 0	
#define WCOREPUSH    wspacepush(ctrl)
#define WCOREPOP     wspacepop(ctrl)
#define UNMATCHED   -1
#define gk_SWAP(a, b, tmp) do {(tmp) = (a); (a) = (b); (b) = (tmp);} while(0) 
#define MAKECSR(i, n, a) \
   do { \
     for (i=1; i<n; i++) a[i] += a[i-1]; \
     for (i=n; i>0; i--) a[i] = a[i-1]; \
     a[0] = 0; \
   } while(0) 
   #define INC_DEC(a, b, val) do {(a) += (val); (b) -= (val);} while(0)

   #define ListInsert(n, lind, lptr, i) \
   do { \
     lind[n] = i; \
     lptr[i] = (n)++;\
   } while(0) 

#define ListDelete(n, lind, lptr, i) \
   do { \
     lind[lptr[i]] = lind[--(n)]; \
     lptr[lind[n]] = lptr[i]; \
     lptr[i] = -1; \
   } while(0) 

#define BNDInsert(nbnd, bndind, bndptr, vtx) \
  ListInsert(nbnd, bndind, bndptr, vtx)

#define BNDDelete(nbnd, bndind, bndptr, vtx) \
  ListDelete(nbnd, bndind, bndptr, vtx)

#define KEY_LT(a, b) ((a) > (b))


typedef struct graph_t {
  int nvtxs, nedges;	
  int ncon;		
  int *xadj;		
  int *vwgt;		
  int *vsize;		
  int *adjncy;        
  int *adjwgt;       
  int *tvwgt;         
  float *invtvwgt;     
  int free_xadj, free_vwgt, free_vsize, free_adjncy, free_adjwgt;
  int *label;
  int *cmap;
  int mincut, minvol;
  int *where, *pwgts;
  int nbnd;
  int *bndptr, *bndind;
  int *id, *ed;
  ckrinfo_t *ckrinfo;   
  vkrinfo_t *vkrinfo;   
  struct graph_t *coarser, *finer;
  //coarsen
  int *cuda_maxvwgt;
  int *cuda_match;
  int *cuda_cmap;
  int *cuda_nvtxs;
  int *cuda_xadj;
  int *cuda_adjncy;
  int *cuda_adjwgt;
  int *cuda_vwgt;
  int *cuda_real_edge;
  int *cuda_real_nvtxs;
  int *cuda_cnvtxs;
  int *cuda_s;
  int *cuda_ccxadjwgt;
  int *cuda_scan_adjwgt_original;
  int *cuda_scan_nedges_original;
  int *cuda_js;
  int *cuda_scan_cnedges_original;
  int *cuda_real_cnedges;
  int *cuda_temp_cnedges;
  //uncoarsen
  int *cuda_maxwgt;
  int *cuda_minwgt;
  int *cuda_where;
  int *cuda_pwgts;
  int *cuda_bnd;
  int *cuda_bndto;
  int *cuda_bndnum;
  int *cpu_bndnum;
  //int *cuda_gain;
  int *bnd_info;
  int *cuda_real_bnd_num;
  int *cuda_real_bnd;
  int *cuda_nparts;
  int *cuda_bn;
  int *cuda_bt;
  int *cuda_g;
  int *cuda_tvwgt;
  float *cuda_ctrl_tpwgts;


} graph_t;

//
int *cu_info;
int *cu_bn;
int *cu_bt;
int *cu_g;
int  *cu_csr;
int *cu_que;
//


typedef struct gk_mop_t {
  int type;
  ssize_t nbytes;
  void *ptr;
} gk_mop_t;
typedef struct gk_mcore_t {
 
  size_t coresize;     
  size_t corecpos;    
  void *core;	       
  size_t nmops;         
  size_t cmop;         
  gk_mop_t *mops;      
  size_t num_callocs;   
  size_t num_hallocs;   
  size_t size_callocs;  
  size_t size_hallocs;  
  size_t cur_callocs;   
  size_t cur_hallocs;  
  size_t max_callocs;   
  size_t max_hallocs;   

} gk_mcore_t;

typedef struct cnbr_t {
  int pid;           
  int ed;          
} cnbr_t;

typedef struct vnbr_t {
  int pid;           
  int ned;            
  int gv;             
} vnbr_t;

typedef struct ctrl_t {
  int CoarsenTo;		
  int nIparts;               
  int no2hop;                 
  int minconn;                
  int contig;                 
  int nseps;		
  int ufactor;                
  int compress;               
  int ccorder;                
  int seed;                   
  int ncuts;                  
  int niter;                 
  int numflag;                
  int *maxvwgt;		
  int ncon;                  
  int nparts;                
  float pfactor;		
  float *ubfactors;            
  float *tpwgts;               
  float *pijbm;               
  float cfactor;               
  double TotalTmr, InitPartTmr, MatchTmr, ContractTmr, CoarsenTmr, UncoarsenTmr, 
         RefTmr, ProjectTmr, SplitTmr, Aux1Tmr, Aux2Tmr, Aux3Tmr;
  gk_mcore_t *mcore;    
  size_t nbrpoolsize;      
  size_t nbrpoolcpos;      
  size_t nbrpoolreallocs;  
  cnbr_t *cnbrpool;    
  vnbr_t *vnbrpool;     
  int *maxnads;               
  int *nads;                  
  int **adids;                
  int **adwgts;               
  int *pvec1, *pvec2;         

} ctrl_t;

typedef struct rkv_t{
  float key;
  int val;
} rkv_t;

typedef struct {
  ssize_t nnodes;
  ssize_t maxnodes;
  rkv_t   *heap;
  ssize_t *locator;
} rpq_t;

int irand() 
{
  if (sizeof(int) <= sizeof(int32_t)) 
    return (int)(uint32_t)rand();
  else 
    return (int)(uint64_t)rand(); 
}

static __thread gk_mcore_t *gkmcore = NULL;

  int  isum(size_t n, int *x)
{

  size_t i;
  int sum = 0;

  for (i=0; i<n; i++, x+=1)
    sum += (*x);

  return sum;
  
}

float rsum(size_t n, float *x)
{
  size_t i;
  float sum = 0;

  for (i=0; i<n; i++, x+=1)
    sum += (*x);

  return sum;
}

float   *rscale(size_t n, float alpha, float *x)
{
  size_t i;

  for (i=0; i<n; i++, x+=1)
    (*x) *= alpha;
  return x;
}



int ComputeCut(graph_t *graph, int *where)
{
  int i, j, cut;


    for (cut=0, i=0; i<graph->nvtxs; i++) {
      for (j=graph->xadj[i]; j<graph->xadj[i+1]; j++)
        if (where[i] != where[graph->adjncy[j]])
          cut += graph->adjwgt[j];
    }
  

  return cut/2;
}


ctrl_t *SetupCtrl(int nparts, float *tpwgts, float *ubvec)
{
  int i;
  ctrl_t *ctrl;
  ctrl = (ctrl_t *)malloc(sizeof(ctrl_t));
  memset((void *)ctrl, 0, sizeof(ctrl_t));
      ctrl->ncuts   = 1;
      ctrl->niter   = 10;
      ctrl->ufactor =1 ;
      ctrl->CoarsenTo = 200;
      ctrl->no2hop  = 0;
      ctrl->seed    = -1;
      ctrl->numflag = 0;
      ctrl->nparts  = nparts; 
      ctrl->maxvwgt = (int*)malloc(sizeof(int));
      ctrl->maxvwgt[0] =0;  
    ctrl->tpwgts = (float*)malloc(sizeof(float)*nparts);
      for (i=0; i<nparts; i++) {
          ctrl->tpwgts[i] = 1.0/nparts;
      }
    ctrl->ubfactors = (float*)malloc(sizeof(float));
    ctrl->ubfactors[0] =1.03;
   ctrl->pijbm =(float*) malloc(sizeof(float)*nparts);
   ctrl->seed=1;
    return ctrl;
    
}


void InitGraph(graph_t *graph) 
{
  memset((void *)graph, 0, sizeof(graph_t));


  graph->nvtxs     = -1;
  graph->nedges    = -1;
  graph->ncon      = -1;
  graph->mincut    = -1;
  graph->minvol    = -1;
  graph->nbnd      = -1;

  graph->xadj      = NULL;
  graph->vwgt      = NULL;
  graph->vsize     = NULL;
  graph->adjncy    = NULL;
  graph->adjwgt    = NULL;
  graph->label     = NULL;
  graph->cmap      = NULL;
  graph->tvwgt     = NULL;
  graph->invtvwgt  = NULL;


  graph->free_xadj   = 1;
  graph->free_vwgt   = 1;
  graph->free_vsize  = 1;
  graph->free_adjncy = 1;
  graph->free_adjwgt = 1;


  graph->where     = NULL;
  graph->pwgts     = NULL;
  graph->id        = NULL;
  graph->ed        = NULL;
  graph->bndptr    = NULL;
  graph->bndind    = NULL;
  graph->ckrinfo   = NULL;
  graph->vkrinfo   = NULL;

  graph->coarser   = NULL;
  graph->finer     = NULL;
}

graph_t *CreateGraph(void)
{
  graph_t *graph;

  graph = (graph_t *)malloc(sizeof(graph_t));

  InitGraph(graph);

  return graph;
}


void SetupGraph_tvwgt(graph_t *graph)
{
 
  if (graph->tvwgt == NULL) 
    graph->tvwgt  = (int*)malloc(sizeof(int));
  if (graph->invtvwgt == NULL) 
    graph->invtvwgt =(float*) malloc(sizeof(float));

    graph->tvwgt[0]    = isum(graph->nvtxs, graph->vwgt);
    
   
    graph->invtvwgt[0] = 1.0/(graph->tvwgt[0] > 0 ? graph->tvwgt[0] : 1);
  
}


void SetupGraph_label(graph_t *graph)
{
  int i;

  if (graph->label == NULL)
    graph->label =(int *) malloc(sizeof(int)*(graph->nvtxs));


  for (i=0; i<graph->nvtxs; i++)
  {
    graph->label[i] = i;
  }
}


graph_t *SetupGraph(ctrl_t *ctrl, int nvtxs, int *xadj, 
             int *adjncy, int *vwgt , int *adjwgt) 
{
  int i;

  graph_t *graph;
  graph = CreateGraph();
  graph->nvtxs  = nvtxs;
  graph->nedges = xadj[nvtxs];
  graph->xadj      = xadj;
  graph->free_xadj = 0;
  graph->adjncy      = adjncy;
  graph->free_adjncy = 0;
  if (vwgt) {
    graph->vwgt      = vwgt;
    graph->free_vwgt = 0;
  }
  else {
    vwgt = graph->vwgt = (int*)malloc(sizeof(int)*nvtxs);
    for(i=0;i<nvtxs;i++)
    {
        vwgt[i]=graph->vwgt[i]=1;
    }
  }
  
     graph->tvwgt    = (int*)malloc(sizeof(int));
    graph->invtvwgt = (float*)malloc(sizeof(float));

    graph->tvwgt[0]    = isum(nvtxs, vwgt);
    

    graph->invtvwgt[0] = 1.0/(graph->tvwgt[0] > 0 ? graph->tvwgt[0] : 1);
    
    if (adjwgt) {
      graph->adjwgt      = adjwgt;
      graph->free_adjwgt = 0;
    }
    else {
      adjwgt = graph->adjwgt=(int*)malloc(sizeof(int)*(graph->nedges));
      for(i=0;i<graph->nedges;i++)
      {
          adjwgt[i] = graph->adjwgt[i]=1;
      }
    }
   SetupGraph_tvwgt(graph);
   SetupGraph_label(graph);


  return graph;
}


gk_mcore_t *gk_mcoreCreate(size_t coresize)
{
  gk_mcore_t *mcore;
  mcore = (gk_mcore_t *)malloc(sizeof(gk_mcore_t));
  memset(mcore, 0, sizeof(gk_mcore_t));
  mcore->coresize = coresize;
  mcore->corecpos = 0;
  mcore->core = (coresize == 0 ? NULL : (size_t*)malloc(sizeof(size_t)*(mcore->coresize)));

  mcore->nmops = 2048;
  mcore->cmop  = 0;
  mcore->mops  = (gk_mop_t *)malloc((mcore->nmops)*sizeof(gk_mop_t));

  return mcore;
}


void AllocateWorkSpace(ctrl_t *ctrl, graph_t *graph)
{
  size_t coresize;


      coresize = 3*(graph->nvtxs+1)*sizeof(int) + 
                 5*(ctrl->nparts+1)*sizeof(int) + 
                 5*(ctrl->nparts+1)*sizeof(float);
  ctrl->mcore = gk_mcoreCreate(coresize);
  ctrl->nbrpoolsize = 0;
  ctrl->nbrpoolcpos = 0;
}

void gk_mcoreAdd(gk_mcore_t *mcore, int type, size_t nbytes, void *ptr)
{
  if (mcore->cmop == mcore->nmops) {
    mcore->nmops *= 2;
    mcore->mops = (gk_mop_t*)realloc(mcore->mops, mcore->nmops*sizeof(gk_mop_t));
    if (mcore->mops == NULL) 
      exit(0);
  }

  mcore->mops[mcore->cmop].type   = type;
  mcore->mops[mcore->cmop].nbytes = nbytes;
  mcore->mops[mcore->cmop].ptr    = ptr;
  mcore->cmop++;

  switch (type) {
    case 1:
      break;

    case 2:
      mcore->num_callocs++;
      mcore->size_callocs += nbytes;
      mcore->cur_callocs  += nbytes;
      if (mcore->max_callocs < mcore->cur_callocs)
        mcore->max_callocs = mcore->cur_callocs;
      break;

    case 3:
      mcore->num_hallocs++;
      mcore->size_hallocs += nbytes;
      mcore->cur_hallocs  += nbytes;
      if (mcore->max_hallocs < mcore->cur_hallocs)
        mcore->max_hallocs = mcore->cur_hallocs;
      break;
    default:
      exit(0);
  }
}


void *gk_mcoreMalloc(gk_mcore_t *mcore, size_t nbytes)
{
  void *ptr;

  nbytes += (nbytes%8 == 0 ? 0 : 8 - nbytes%8);

  if (mcore->corecpos + nbytes < mcore->coresize) {

    ptr = ((char *)mcore->core)+mcore->corecpos;
    mcore->corecpos += nbytes;
    gk_mcoreAdd(mcore, 2, nbytes, ptr);
  }
  else {

    ptr = (size_t*)malloc(nbytes);

    gk_mcoreAdd(mcore, 3, nbytes, ptr);
  }
  return ptr;
}


void *wspacemalloc(ctrl_t *ctrl, size_t nbytes)
{
  return gk_mcoreMalloc(ctrl->mcore, nbytes);
}


float *rwspacemalloc(ctrl_t *ctrl)
{
  return (float *)wspacemalloc(ctrl, 2*sizeof(float));
}


void Setup2WayBalMultipliers(ctrl_t *ctrl, graph_t *graph, float *tpwgts)
{
  int i;

  for (i=0; i<2; i++) {
      ctrl->pijbm[i] = graph->invtvwgt[0]/tpwgts[i];
  }
}


int  *icopy(size_t n, int *a, int *b)
{
  return (int *)memmove((void *)b, (void *)a, sizeof(int)*n);
}


int *iwspacemalloc(ctrl_t *ctrl, size_t n)
{
  return (int *)wspacemalloc(ctrl, n*sizeof(int));
}


int *iset(size_t n, int val, int *x)
{
  size_t i;
  for (i=0; i<n; i++)
    x[i] = val;
  return x;
}


int irandInRange(int max) 
{
  return (int)((irand())%max); 
}


void irandArrayPermute(int n, int *p, int nshuffles, int flag)
{
  int i, u, v;
  int tmp;

  if (flag == 1) {

    for (i=0; i<n; i++)
      p[i] = (int)i;
  }

  if (n < 10) {

    for (i=0; i<n; i++) {
      v = irandInRange(n);
      u = irandInRange(n);
      gk_SWAP(p[v], p[u], tmp);
    }
  }
  else {

    for (i=0; i<nshuffles; i++) {
      v = irandInRange(n-3);
      u = irandInRange(n-3);
      gk_SWAP(p[v+0], p[u+2], tmp);
      gk_SWAP(p[v+1], p[u+3], tmp);
      gk_SWAP(p[v+2], p[u+0], tmp);
      gk_SWAP(p[v+3], p[u+1], tmp);
    }
  }
}

void BucketSortKeysInc(ctrl_t *ctrl, int n, int max, int *keys, 
         int *tperm, int *perm)
{
  int i, ii;
  int *counts;
  counts = iset(max+2, 0, iwspacemalloc(ctrl, max+2));
  for (i=0; i<n; i++)
    counts[keys[i]]++;
  MAKECSR(i, max+1, counts);

  for (ii=0; ii<n; ii++) {
    i = tperm[ii];
    perm[counts[keys[i]]++] = i;
  }
}


graph_t *SetupCoarseGraph(graph_t *graph, int cnvtxs)
{
  graph_t *cgraph;

  cgraph = CreateGraph();

  cgraph->nvtxs = cnvtxs;

  cgraph->finer  = graph;
  graph->coarser = cgraph;

  cgraph->xadj     = (int*)malloc(sizeof(int)*(cnvtxs+1));
  cgraph->adjncy   = (int*)malloc(sizeof(int)*(graph->nedges));
  cgraph->adjwgt   = (int*)malloc(sizeof(int)*(graph->nedges));
  cgraph->vwgt     = (int*)malloc(sizeof(int)*cnvtxs);
  cgraph->tvwgt    = (int*)malloc(sizeof(int));
  cgraph->invtvwgt =(float*)malloc(sizeof(float)); 

  return cgraph;
}

graph_t *SetupCoarseGraph1(graph_t *graph, int cnvtxs)
{
  graph_t *cgraph;

  cgraph = CreateGraph();

  cgraph->nvtxs = cnvtxs;

  cgraph->finer  = graph;
  graph->coarser = cgraph;

  cgraph->xadj     = (int*)malloc(sizeof(int)*(cnvtxs+1));

 
  cgraph->tvwgt    = (int*)malloc(sizeof(int));
  cgraph->invtvwgt =(float*)malloc(sizeof(float)); 


  return cgraph;
}




void CreateCoarseGraphNoMask0(ctrl_t *ctrl, graph_t *graph, int cnvtxs, 
         int *match,int level)
{

  int i,j, k, m, istart, iend, nvtxs, nedges, cnedges, v, u;
  int *xadj, *vwgt, *adjncy, *adjwgt;
  int *cmap, *htable;
  int *cxadj, *cvwgt, *cadjncy, *cadjwgt;
  graph_t *cgraph;
 
  nvtxs   = graph->nvtxs;
  xadj    = graph->xadj;
  vwgt    = graph->vwgt;
  adjncy  = graph->adjncy;
  adjwgt  = graph->adjwgt;
  cmap    = graph->cmap;                  
  cgraph = SetupCoarseGraph(graph, cnvtxs);            
  cxadj    = cgraph->xadj;
  cvwgt    = cgraph->vwgt;

  cadjncy  = cgraph->adjncy;
  cadjwgt  = cgraph->adjwgt;                               
  htable = iset(cnvtxs, -1, iwspacemalloc(ctrl, cnvtxs));      

  cxadj[0] = cnvtxs = cnedges = 0; 
  nedges=graph->nedges;
  for (v=0; v<nvtxs; v++) 
  {

    if ((u = match[v]) < v)         
      continue;               

      cvwgt[cnvtxs] = vwgt[v];                 
    nedges = 0;                                                    
    istart = xadj[v];
    iend   = xadj[v+1];                   

    for (j=istart; j<iend; j++)             
    {
         k = cmap[adjncy[j]];                         
         if ((m = htable[k]) == -1)                     
         {
           cadjncy[nedges] = k;    
                               
           cadjwgt[nedges] = adjwgt[j];                      
           htable[k] = nedges++;                                    
         }
         else 
         {
           cadjwgt[m] += adjwgt[j];                                 
         }
    }
   
    if (v != u)                      
     { 
   
           cvwgt[cnvtxs] += vwgt[u];                   
         istart = xadj[u];                                    
         iend   = xadj[u+1];                           
         for (j=istart; j<iend; j++) 
         {
             k = cmap[adjncy[j]];
             if ((m = htable[k]) == -1) 
              {
               cadjncy[nedges] = k;
               cadjwgt[nedges] = adjwgt[j];
               htable[k] = nedges++;
              }
             else 
             {
               cadjwgt[m] += adjwgt[j];
             }
      }

         if ((j = htable[cnvtxs]) != -1)                            
         {
           cadjncy[j]        = cadjncy[--nedges];
           cadjwgt[j]        = cadjwgt[nedges];
           htable[cnvtxs] = -1;
         }
    }


     for (j=0; j<nedges; j++)
       htable[cadjncy[j]] = -1;  

     cnedges         += nedges;
     cxadj[++cnvtxs]  = cnedges;
     cadjncy         += nedges;                                                                 
     cadjwgt         += nedges;

  }


  cgraph->nedges = cnedges;


    cgraph->tvwgt[0]    = isum(cgraph->nvtxs, cgraph->vwgt); 
    
                                      
    cgraph->invtvwgt[0] = 1.0/(cgraph->tvwgt[0] > 0 ? cgraph->tvwgt[0] : 1);           
    


}



int Match_SHEM0(ctrl_t *ctrl, graph_t *graph,int level)
{

  int i, j, pi,k, nvtxs,  cnvtxs, maxidx, maxwgt, 
     last_unmatched,avgdegree,u;
  int *xadj, *vwgt, *adjncy, *adjwgt, *maxvwgt;
  int *match, *cmap, *degrees, *perm, *tperm;

  
  nvtxs  = graph->nvtxs;
  xadj   = graph->xadj;
  vwgt   = graph->vwgt;
  adjncy = graph->adjncy;
  adjwgt = graph->adjwgt;
  cmap   = graph->cmap;
  maxvwgt  = ctrl->maxvwgt;
  
  match   = iset(nvtxs, UNMATCHED, iwspacemalloc(ctrl, nvtxs));
  
  perm    = iwspacemalloc(ctrl, nvtxs);
  tperm   = iwspacemalloc(ctrl, nvtxs);
  degrees = iwspacemalloc(ctrl, nvtxs);         
  
  irandArrayPermute(nvtxs, tperm, nvtxs/8, 1);   

  avgdegree = 0.7*(xadj[nvtxs]/nvtxs);

  for (i=0; i<nvtxs; i++) 
    degrees[i] = (xadj[i+1]-xadj[i] > avgdegree ? avgdegree : xadj[i+1]-xadj[i]);


  BucketSortKeysInc(ctrl, nvtxs, avgdegree, degrees, tperm, perm);         

  
  cnvtxs=0;
  last_unmatched=0;

  for (pi=0; pi<nvtxs; pi++) 
  {
    i = perm[pi];  
    if (match[i] == UNMATCHED) 
    {  

      maxidx = i;                                                                               
      maxwgt = -1;                                                                      

            for (j=xadj[i]; j<xadj[i+1]; j++) 
            {
              k = adjncy[j];
              if (match[k] == UNMATCHED && maxwgt < adjwgt[j] && vwgt[i]+vwgt[k] <= maxvwgt[0]) 
                  {
                       maxidx = k;
                        maxwgt = adjwgt[j];
              }
   
            if (maxidx == i && 3*vwgt[i] < maxvwgt[0]) 
            { 
              maxidx = UNMATCHED;
            }
          }

      
  if (maxidx != UNMATCHED) 
      {
       cmap[i]  = cmap[maxidx] = cnvtxs++;              
        match[i] = maxidx;                                        
        match[maxidx] = i; 
      }
  }
 
  }


                  
  for (cnvtxs=0, i=0; i<nvtxs; i++) 
  {

    if (match[i] == UNMATCHED) {
      match[i] = i;
      cmap[i] = cnvtxs++;                                                    
    }
    else {
      if (i <= match[i]) 
        cmap[i] = cmap[match[i]] = cnvtxs++;
    }

  }


  CreateCoarseGraphNoMask0(ctrl, graph, cnvtxs, match,level);

  return cnvtxs;
}



void Mallocoriginal_coarseninfo(ctrl_t *ctrl,graph_t *graph)
{
    int nvtxs=graph->nvtxs;
    int nedges=graph->nedges;

    cudaMalloc((void**)&graph->cuda_nvtxs, sizeof(int));
    cudaMemcpy(graph->cuda_nvtxs,&graph->nvtxs,sizeof(int),cudaMemcpyHostToDevice);

    cudaMalloc((void**)&graph->cuda_match, nvtxs*sizeof(int));

    cudaMalloc((void**)&graph->cuda_xadj, (nvtxs+1)*sizeof(int));
    cudaMemcpy(graph->cuda_xadj,graph->xadj,(nvtxs+1)*sizeof(int),cudaMemcpyHostToDevice);
   
    cudaMalloc((void**)&graph->cuda_vwgt,  (nvtxs+1)*sizeof(int));
    cudaMemcpy(graph->cuda_vwgt,graph->vwgt,nvtxs*sizeof(int),cudaMemcpyHostToDevice);
    
    cudaMalloc((void**)&graph->cuda_adjncy,  nedges*sizeof(int));
    cudaMemcpy(graph->cuda_adjncy,graph->adjncy,nedges*sizeof(int),cudaMemcpyHostToDevice);
    
    cudaMalloc((void**)&graph->cuda_adjwgt,  nedges*sizeof(int));
    cudaMemcpy(graph->cuda_adjwgt,graph->adjwgt,nedges*sizeof(int),cudaMemcpyHostToDevice);

    cudaMalloc((void**)&graph->cuda_cmap, nvtxs*sizeof(int));

    cudaMalloc((void**)&graph->cuda_maxvwgt, sizeof(int));
    
     cudaMalloc((void**)&graph->cuda_s, (nvtxs)* sizeof(int));

     cudaMalloc((void**)&graph->cuda_cnvtxs, sizeof(int));

     cudaMalloc((void**)&graph->cuda_scan_nedges_original,  (graph->nedges)*sizeof(int));

     cudaMalloc((void**)&graph->cuda_scan_cnedges_original,  (graph->nedges)*sizeof(int));

     cudaMalloc((void**)&graph->cuda_scan_adjwgt_original,  (graph->nedges)*sizeof(int));

     cudaMalloc((void**)&graph->cuda_ccxadjwgt, (graph->nedges)* sizeof(int));

     cudaMalloc((void**)&graph->cuda_temp_cnedges, (graph->nedges)* sizeof(int));




}


  void Malloccoarseninfo(ctrl_t *ctrl,graph_t *graph)
{
    int nvtxs=graph->nvtxs;
    int nedges=graph->nedges;

    cudaMalloc((void**)&graph->cuda_nvtxs, sizeof(int));
    cudaMemcpy(graph->cuda_nvtxs,&graph->nvtxs,sizeof(int),cudaMemcpyHostToDevice);

    cudaMalloc((void**)&graph->cuda_match, nvtxs*sizeof(int));


    cudaMalloc((void**)&graph->cuda_cmap, nvtxs*sizeof(int));

    cudaMalloc((void**)&graph->cuda_maxvwgt, sizeof(int));
    
     cudaMalloc((void**)&graph->cuda_s, (nvtxs)* sizeof(int));

     cudaMalloc((void**)&graph->cuda_cnvtxs, sizeof(int));

     cudaMalloc((void**)&graph->cuda_scan_nedges_original,  (graph->nedges)*sizeof(int));

     cudaMalloc((void**)&graph->cuda_scan_cnedges_original,  (graph->nedges)*sizeof(int));

     cudaMalloc((void**)&graph->cuda_scan_adjwgt_original,  (graph->nedges)*sizeof(int));

     cudaMalloc((void**)&graph->cuda_ccxadjwgt, (graph->nedges)* sizeof(int));

     cudaMalloc((void**)&graph->cuda_temp_cnedges, (graph->nedges)* sizeof(int));
  }





__global__ void initcuda_match(int *cuda_match)
{
  int ii;
  ii=blockIdx.x*blockDim.x+threadIdx.x;
  cuda_match[ii]=-1;
}


__global__ void cuda_shem(  int *cuda_nvtxs,  int *match, int *xadj,int *vwgt,  int *adjwgt,  \
  int *adjncy,  int *maxvwgt)
{
  
  int  pi;
  int ii,iii;
  int i,j,k,maxidx,maxwgt,u;
  ii=blockIdx.x*blockDim.x+threadIdx.x,iii= blockIdx.x*blockDim.x+threadIdx.x+1;
  int b_start,b_end;
  int tt=1024;
  int nvtxs=cuda_nvtxs[0];
  if(nvtxs%tt==0)
  {
    b_start=ii*(nvtxs/tt);
    b_end=(ii+1)*(nvtxs/tt);
  }
  else
  {
    int b=nvtxs/tt;
    int a=b+1;
    int x=nvtxs-b*tt;
    int y=tt-x;
    if(ii<x)
    {
      b_start=ii*a;
      b_end=(ii+1)*a;
    }
    else
    {
      b_start=x*a+(ii-x)*b;
      b_end=x*a+(ii+1-x)*b;
    }
  }


  
  for (pi=b_start; pi<b_end; pi++) 
  {
     i=pi;
    if (match[i] == -1) 
    {  

      maxidx = i;                                                                               
      maxwgt = -1;                                                                      

            for (j=xadj[i]; j<xadj[i+1]; j++) 
            {
              k = adjncy[j];
              if (match[k] == -1 && maxwgt < adjwgt[j] && vwgt[i]+vwgt[k] <= maxvwgt[0]) 
                  {
                       maxidx = k;
                        maxwgt = adjwgt[j];
              }
   
            if (maxidx == i && 3*vwgt[i] < maxvwgt[0]) 
            { 
              maxidx = -1;
            }
          }
  if (maxidx != -1) 
      {    
        match[i] = maxidx;  
         atomicExch(&match[maxidx],i);                                 
      }
     }


  
     }
   
}


__global__ void cuda_shem1(  int *cuda_nvtxs,  int *match, int *xadj,int *vwgt,  int *adjwgt,  \
  int *adjncy,  int *maxvwgt)
{
  
  int  pi;
  int ii,iii;
  int i,j,k,maxidx,maxwgt,u;
  ii=blockIdx.x*blockDim.x+threadIdx.x,iii= blockIdx.x*blockDim.x+threadIdx.x+1;
  int b_start,b_end;
   int nvtxs=cuda_nvtxs[0];
  
   for (pi=0; pi<nvtxs; pi++) 
  {
    i = pi;  
    if (match[i] == -1) 
    {  

      maxidx = i;                                                                               
      maxwgt = -1;                                                                      

            for (j=xadj[i]; j<xadj[i+1]; j++) 
            {
        
              if (match[k] == -1 && maxwgt < adjwgt[j] && vwgt[i]+vwgt[k] <= maxvwgt[0]) 
                  {
                       maxidx = k;
                        maxwgt = adjwgt[j];
              }
    

            if (maxidx == i && 3*vwgt[i] < maxvwgt[0]) 
            { 
              maxidx = -1;
            }
         
          }

      
  if (maxidx != -1) 
      {            
        match[i] = maxidx;  
        atomicExch(&match[maxidx],i);   
      }
  }
 
  }
   
}


__global__ void cuda_cleanv(int *match,int *s)
{
  int ii,iii,pi,u;
  pi=blockIdx.x*blockDim.x+threadIdx.x;
  

  
        s[pi]=1;
         if(match[pi]!=-1)
         {
           u=match[pi];
           if(match[u]!=pi)
             s[pi]=0;
         }

  
} 


__global__ void findc(int *cuda_match,int *cuda_cmap,int *cuda_cnvtxs,int *s)
{

  int  pi;
  int ii,iii;
  int i,j;
  
  pi=blockIdx.x*blockDim.x+threadIdx.x;

  
 if(s[pi]==0)
 cuda_match[pi]=pi;

 
 if (cuda_match[pi] == -1) 
  {

    cuda_match[pi] = pi;  
    cuda_cmap[pi] = atomicAdd(&cuda_cnvtxs[0],1);  
                     
  }
  else  if(pi<=cuda_match[pi])
  {
      cuda_cmap[pi] =cuda_cmap[cuda_match[pi]]= atomicAdd(&cuda_cnvtxs[0],1); 

  }


}


__global__ void findscaninfo(int *match,int *real,int *cmap,int *cnedges,int *xadj)
{
  int ii,iii,pi,u;
  pi=blockIdx.x*blockDim.x+threadIdx.x;


    if ((u = match[pi]) >= pi)           
    {
      real[cmap[pi]]=pi;

      if(u!=pi)
      cnedges[cmap[pi]]=(xadj[pi+1]-xadj[pi])+(xadj[u+1]-xadj[u]);
      else
      cnedges[cmap[pi]]=(xadj[pi+1]-xadj[pi]);


    }
  
}

__global__ void initcudajs(int *cuda_js)
{
  int ii;
  ii=blockIdx.x*blockDim.x+threadIdx.x;
  cuda_js[ii]=0;
}


__global__ void find_cnvtxsedge_original(int *cuda_scan_nedges,int *cuda_scan_order,int *cuda_xadj,\
  int *cuda_match,int *cuda_adjncy,int *cuda_scan_nedges_original,int *cuda_cmap,int *cvwgt,int *vwgt,int *js,\
  int *cuda_scan_cnedges_original,int *cuda_scan_adjwgt_original,int *cuda_adjwgt,int *cuda_ccxadjwgts)
  {
    int ii,iii,pi,u,v,j,istart,iend,k,nedges,m,jjs,i;
    
    pi=blockIdx.x*blockDim.x+threadIdx.x;

    int pp,ppp;
    int pii;
   
    pii=cuda_scan_order[pi];
    
    pp=cuda_scan_nedges[cuda_cmap[pii]];
    ppp=cuda_scan_nedges[cuda_cmap[pii]+1];

    u = cuda_match[pii];
 
      for(i=cuda_xadj[pii];i<cuda_xadj[pii+1];i++)
      {

        cuda_scan_nedges_original[pp+js[cuda_cmap[pii]]]=cuda_adjncy[i];
        cuda_scan_cnedges_original[pp+js[cuda_cmap[pii]]]=cuda_cmap[cuda_adjncy[i]];
        cuda_scan_adjwgt_original[pp+js[cuda_cmap[pii]]]=cuda_adjwgt[i];
        cuda_ccxadjwgts[pp+js[cuda_cmap[pii]]]=0;
        js[cuda_cmap[pii]]++;
      } 
  
  
      if(u!=pii)
      {
      for(i=cuda_xadj[u];i<cuda_xadj[u+1];i++)
      {
        cuda_scan_nedges_original[pp+js[cuda_cmap[pii]]]=cuda_adjncy[i];
        cuda_scan_cnedges_original[pp+js[cuda_cmap[pii]]]=cuda_cmap[cuda_adjncy[i]];
        cuda_scan_adjwgt_original[pp+js[cuda_cmap[pii]]]=cuda_adjwgt[i];
        cuda_ccxadjwgts[pp+js[cuda_cmap[pii]]]=0;
          js[cuda_cmap[pii]]++;
       
      }     
      }
   
      if(u!=pii)
    {
      cvwgt[cuda_cmap[pii]]=vwgt[pii]+vwgt[u];
    }
    else
    {
      cvwgt[cuda_cmap[pii]]=vwgt[pii];
    }



   
    

}

__global__ void find_cnvtxsedge_original1(int *cuda_scan_nedges,int *cuda_scan_order,int *cuda_xadj,\
  int *cuda_match,int *cuda_adjncy,int *cuda_scan_nedges_original,int *cuda_cmap,int *cvwgt,int *vwgt,int *js,\
  int *cuda_scan_cnedges_original,int *cuda_scan_adjwgt_original,int *cuda_adjwgt,int *cuda_ccxadjwgts)
  {
    int ii,iii,pi,u,v,j,istart,iend,k,nedges,m,jjs,i;
    
    pi=blockIdx.x*blockDim.x+threadIdx.x;

    int pp,ppp;
    int pii;
    u=cuda_match[pi];
    
    if(pi>u)
    {
      pp=cuda_scan_nedges[cuda_cmap[pi]]+cuda_xadj[u+1]-cuda_xadj[u];
    }
    else
    {
      pp=cuda_scan_nedges[cuda_cmap[pi]];
    }

    jjs=0;
    for(i=cuda_xadj[pi];i<cuda_xadj[pi+1];i++)
    {
        int pt=pp+jjs;
      cuda_scan_nedges_original[pt]=cuda_adjncy[i];
      cuda_scan_cnedges_original[pt]=cuda_cmap[cuda_adjncy[i]];
      cuda_scan_adjwgt_original[pt]=cuda_adjwgt[i];
      cuda_ccxadjwgts[pt]=0;
      jjs++;
    } 

    if(u!=pi)
    {
      cvwgt[cuda_cmap[pi]]=vwgt[pi]+vwgt[u];
    }
    else
    {
      cvwgt[cuda_cmap[pi]]=vwgt[pi];
    }


}



__device__ void quicksort(int *arr,int low ,int high,int *arr2)
{
   if(low<high)
   {
      int i=low;
      int j=high;
      int k=arr[low];
      int k2=arr2[low];
      while(i<j)
      {
          while(i<j&&arr[j]>=k)
          {
            j--;
          }
          if(i<j)
          {
            arr[i]=arr[j];
            arr2[i]=arr2[j];
            i++;
          }
          while(i<j&&arr[i]<k)
          {
            i++;
          }
          if(i<j)
          {
            arr[j]=arr[i];
            arr2[j]=arr2[i];
            j--;
          }
      }
      arr[i]=k;
      arr2[i]=k2;
      quicksort(arr,low,i-1,arr2);
      quicksort(arr,i+1,high,arr2);
   }
}

__device__ void swap(int *x,int *y)
{
  int t=*x;
  *x=*y;
  *y=t;
}

__device__ void quick_sort(int *arr,int start,int end,int *arr2)
{
  if(start>=end)
  return;

  int mid=arr[end];
  int left=start,right=end-1;
  while(left<right)
  {
    while(arr[left]<mid&&left<right)
    left++;
    while(arr[right]>=mid&&left<right)
    right--;

    swap(&arr[left],&arr[right]);
    swap(&arr2[left],&arr2[right]);
  }
  if(arr[left]>=arr[end])
  {
    swap(&arr[left],&arr[end]);
    swap(&arr2[left],&arr2[end]);
  }
  else left++;
  if(left)
    quick_sort(arr,start,left-1,arr2);

    quick_sort(arr,left+1,end,arr2);
}

__device__ void bubble_sort(int *arr,int n1,int n2,int *arr2)
{
  int i,j,temp,temp1;
  for(i=n1;i<n2-1;i++)
     for(j=n1;j<n2-1-i;j++) 
       if(arr[j]>arr[j+1])
       {
         temp=arr[j];
         temp1=arr2[j];
         arr[j]=arr[j+1];
         arr2[j]=arr2[j+1];
         arr[j+1]=temp;
         arr2[j+1]=temp1;
       }
}

__global__ void  Sort_cnedges1(int *cuda_scan_cnedges_original,int *cuda_real_cnedges,\
  int *cuda_scan_nedges,int *cuda_scan_order,int *cuda_cmap,int *cuda_temp_cnedges,int *cuda_scan_adjwgt_original,\
  int *cuda_ccxadjwgt)
  {
      int ii,iii,pi,u,v,j,istart,iend,k,nedges,m,jjs,i;
      pi=blockIdx.x*blockDim.x+threadIdx.x;
  
      int pp,ppp;
      int pii;
    
      
      pii=cuda_scan_order[pi];
      
      pp=cuda_scan_nedges[cuda_cmap[pii]];
      ppp=cuda_scan_nedges[cuda_cmap[pii]+1];
 

      quicksort(cuda_scan_cnedges_original,pp,ppp-1,cuda_scan_adjwgt_original);
     
  }




__global__ void  Sort_cnedges22(int *cuda_scan_cnedges_original,int *cuda_real_cnedges,\
  int *cuda_scan_nedges,int *cuda_scan_order,int *cuda_cmap,int *cuda_temp_cnedges,int *cuda_scan_adjwgt_original,\
  int *cuda_ccxadjwgt)
{
    
      int ii,iii,pi,u,v,j,istart,iend,k,nedges,m,jss,i;
      pi=blockIdx.x*blockDim.x+threadIdx.x;
  
      int pp,ppp;
      int pii;
  
      pii=cuda_scan_order[pi];
      pp=cuda_scan_nedges[cuda_cmap[pii]];
      ppp=cuda_scan_nedges[cuda_cmap[pii]+1];
    
    jss=0;
    for(i=pp;i<ppp;i++)
    {
      if(cuda_scan_cnedges_original[i]==cuda_cmap[pii])
      {
      continue;
      }
      else
      {
           if((cuda_scan_cnedges_original[i]!=cuda_scan_cnedges_original[i-1])) //越界了
            {
              cuda_temp_cnedges[pp+jss]=cuda_scan_cnedges_original[i];
              cuda_ccxadjwgt[pp+jss]=cuda_scan_adjwgt_original[i];
              jss++;
            }
            else
            {
              cuda_ccxadjwgt[pp+jss-1]+=cuda_scan_adjwgt_original[i];
            }

      }
    }

    cuda_real_cnedges[cuda_cmap[pii]]=jss;

      
}

__global__ void  Sort_cnedges2(int *cuda_scan_cnedges_original,int *cuda_real_cnedges,\
  int *cuda_scan_nedges,int *cuda_scan_order,int *cuda_cmap,int *cuda_temp_cnedges,int *cuda_scan_adjwgt_original,\
  int *cuda_ccxadjwgt)
{
    
      int ii,iii,pi,u,v,j,istart,iend,k,nedges,m,jss,i;
      pi=blockIdx.x*blockDim.x+threadIdx.x;
  
      int pp,ppp;
      int pii;
  
      pii=cuda_scan_order[pi];

      jss=0;
      
      pp=cuda_scan_nedges[cuda_cmap[pii]];
      ppp=cuda_scan_nedges[cuda_cmap[pii]+1];
  
   
     if(pp!=ppp)
     {
       if(pp==ppp-1)
      {
            if(cuda_scan_cnedges_original[pp]!=cuda_cmap[pii])
            {
               cuda_temp_cnedges[pp+jss]=cuda_scan_cnedges_original[pp];
               cuda_ccxadjwgt[pp+jss]+=cuda_scan_adjwgt_original[pp];
               jss++;
            }
      }
     
      else
      { 
         if(cuda_scan_cnedges_original[pp]!=cuda_cmap[pii])
        {
          cuda_temp_cnedges[pp+jss]=cuda_scan_cnedges_original[pp];
          cuda_ccxadjwgt[pp+jss]+=cuda_scan_adjwgt_original[pp];
          jss++;
         }
         
  
         for(i=pp+1;i<ppp;i++)
          {
            if(cuda_scan_cnedges_original[i]==cuda_cmap[pii])
            continue;
            else if((cuda_scan_cnedges_original[i]==cuda_scan_cnedges_original[i-1]))
            {
              cuda_ccxadjwgt[pp+jss-1]+=cuda_scan_adjwgt_original[i];
            }
            else
            {
              cuda_temp_cnedges[pp+jss]=cuda_scan_cnedges_original[i];
              cuda_ccxadjwgt[pp+jss]=cuda_scan_adjwgt_original[i];
              jss++;
            }
          }
  
      }

      cuda_real_cnedges[cuda_cmap[pii]]=jss;
  
    }
  
}




__global__ void Create_cadjncy(int *cuda_scan_nedges,int *cuda_scan_order, int *cuda_adjwgt,\
  int *cuda_cadjncy,int *cuda_cadjwgt,int *cuda_cmap,int *cuda_real_cnedges,int *cuda_cxadj,int *cuda_temp_cnedges,int *cuda_ccxadjwgt)
  {
    int ii,iii,pi,u,v,j,istart,iend,k,nedges,m,jjs,i;
    pi=blockIdx.x*blockDim.x+threadIdx.x;

    int pp,ppp;
    int pii;
    
  
    pii=cuda_scan_order[pi];
    
    pp=cuda_scan_nedges[cuda_cmap[pii]];
    
    for(i=0;i<cuda_real_cnedges[cuda_cmap[pii]];i++)
    {
   
       cuda_cadjncy[cuda_cxadj[cuda_cmap[pii]]+i]=cuda_temp_cnedges[pp+i];
       cuda_cadjwgt[cuda_cxadj[cuda_cmap[pii]]+i]=cuda_ccxadjwgt[pp+i];

    }
    

     



  }


void Freecudacoarsen(graph_t *graph)
{

  cudaFree(graph->cuda_maxvwgt);
  cudaFree(graph->cuda_match);
  cudaFree(graph->cuda_real_edge);
  cudaFree(graph->cuda_real_nvtxs);
  cudaFree(graph->cuda_cnvtxs);
  cudaFree(graph->cuda_s);
  cudaFree(graph->cuda_ccxadjwgt);
  cudaFree(graph->cuda_scan_adjwgt_original);
  cudaFree(graph->cuda_scan_nedges_original);
  cudaFree(graph->cuda_js);
  cudaFree(graph->cuda_scan_cnedges_original);

  cudaFree(graph->cuda_temp_cnedges);

}



void CreateCoarseGraph(ctrl_t *ctrl, graph_t *graph, int cnvtxs, int level,int *scan_edge)
{
    int nvtxs=graph->nvtxs;
    int nedges=graph->nedges;
    graph_t *cgraph;

    cgraph = SetupCoarseGraph1(graph, cnvtxs);  
    
    int length=cnvtxs+1;
   
   t_cck1_start=clock(); 
   cudaDeviceSynchronize();
    thrust::exclusive_scan(scan_edge,scan_edge+length,scan_edge);
    cudaDeviceSynchronize();
      t_cck1_end=clock(); 
    duration_cck1 += (double)(t_cck1_end - t_cck1_start) / CLOCKS_PER_SEC/test_time; 
 
   

    cudaMemcpy(  graph->cuda_real_edge, scan_edge, (cnvtxs+1)* sizeof(int), cudaMemcpyHostToDevice);

    cudaMalloc((void**)&cgraph->cuda_vwgt, cnvtxs*sizeof(int));
    
    cudaMalloc((void**)&graph->cuda_js, cnvtxs*sizeof(int));

    cudaMalloc((void**)&graph->cuda_real_cnedges, (cnvtxs+1)* sizeof(int));
    cudaDeviceSynchronize();
    initcudajs<<<cnvtxs,1>>>(graph->cuda_js);
    cudaDeviceSynchronize();

  
    t_cck2_start=clock(); 
    cudaDeviceSynchronize();


     find_cnvtxsedge_original1<<<nvtxs,1>>>(graph->cuda_real_edge,graph->cuda_real_nvtxs,graph->cuda_xadj,graph->cuda_match,\
        graph->cuda_adjncy,graph->cuda_scan_nedges_original,graph->cuda_cmap,cgraph->cuda_vwgt,graph->cuda_vwgt,graph->cuda_js,\
        graph->cuda_scan_cnedges_original,graph->cuda_scan_adjwgt_original,graph->cuda_adjwgt,graph->cuda_ccxadjwgt);


      cudaDeviceSynchronize();
      t_cck2_end=clock(); 
      duration_cck2 += (double)(t_cck2_end - t_cck2_start) / CLOCKS_PER_SEC/test_time; 

   
        t_cck3_start=clock(); 
    cudaDeviceSynchronize();

   bb_segsort(graph->cuda_scan_cnedges_original, graph->cuda_scan_adjwgt_original, graph->nedges, graph->cuda_real_edge, cnvtxs);
   cudaDeviceSynchronize();
   t_cck3_end=clock(); 
   duration_cck3 += (double)(t_cck3_end - t_cck3_start) / CLOCKS_PER_SEC/test_time; 



   t_cck4_start=clock(); 
     cudaDeviceSynchronize();
     Sort_cnedges2<<<cnvtxs,1>>>(graph->cuda_scan_cnedges_original,graph->cuda_real_cnedges,graph->cuda_real_edge,graph->cuda_real_nvtxs,\
      graph->cuda_cmap,graph->cuda_temp_cnedges,graph->cuda_scan_adjwgt_original,graph->cuda_ccxadjwgt);
      cudaDeviceSynchronize();
      t_cck4_end=clock(); 
      duration_cck4 += (double)(t_cck4_end - t_cck4_start) / CLOCKS_PER_SEC/test_time; 

    int *cxadj=(int *)malloc(sizeof(int)*(cnvtxs+1));
    cudaMemcpy(  cxadj, graph->cuda_real_cnedges, (cnvtxs+1)* sizeof(int), cudaMemcpyDeviceToHost); 
    
    t_cck5_start=clock(); 
    cudaDeviceSynchronize();
    thrust::exclusive_scan(cxadj,cxadj+length,cxadj);
    cudaDeviceSynchronize();
    t_cck5_end=clock(); 
    duration_cck5 += (double)(t_cck5_end - t_cck5_start) / CLOCKS_PER_SEC/test_time; 

    cgraph->nvtxs=cnvtxs;
    cgraph->nedges=cxadj[cnvtxs];


    cudaMalloc((void**)&cgraph->cuda_xadj, (cnvtxs+1)*sizeof(int));
    cudaMemcpy(  cgraph->cuda_xadj,cxadj,  (cnvtxs+1)* sizeof(int), cudaMemcpyHostToDevice);

    cudaMalloc((void**)&cgraph->cuda_adjncy,   cgraph->nedges*sizeof(int));
      cudaMalloc((void**)&cgraph->cuda_adjwgt,   cgraph->nedges*sizeof(int));

      t_cck6_start=clock(); 
      cudaDeviceSynchronize();
    Create_cadjncy<<<cnvtxs,1>>>(graph->cuda_real_edge,graph->cuda_real_nvtxs,graph->cuda_adjwgt,cgraph->cuda_adjncy,\
      cgraph->cuda_adjwgt,graph->cuda_cmap,graph->cuda_real_cnedges,cgraph->cuda_xadj,\
      graph->cuda_temp_cnedges,graph->cuda_ccxadjwgt);
      cudaDeviceSynchronize();
      t_cck6_end=clock(); 
      duration_cck6 += (double)(t_cck6_end - t_cck6_start) / CLOCKS_PER_SEC/test_time; 

      cgraph->tvwgt[0]    =graph->tvwgt[0];     
      if(level!=0)
      Freecudacoarsen(graph);

}

__global__ void findc1(int *cuda_match,int *cuda_cmap,int *cuda_cnvtxs,int *s)
{

  int  pi;
  int ii,iii;
  int i,j;
  
  pi=blockIdx.x*blockDim.x+threadIdx.x;

  
 if(s[pi]==0||cuda_match[pi] == -1)
 cuda_match[pi]=pi;


}


__global__ void findc2(int *cuda_match,int *cuda_cmap,int *cuda_cnvtxs,int *s)
{

  int  pi;
  int ii,iii;
  int i,j;
  
  pi=blockIdx.x*blockDim.x+threadIdx.x;

  
 if(pi<=cuda_match[pi])
 cuda_cmap[pi]=1;

else
cuda_cmap[pi]=0;


}

__global__ void findc2_5(int *cuda_temp,int *cuda_cmap,int *nvtxs)
{

  int i=nvtxs[0]-1;
  cuda_temp[0]=cuda_cmap[i];
}

__global__ void findc3(int *cuda_match,int *cuda_cmap,int *cuda_cnvtxs,int *s)
{

  int  pi;
  int ii,iii;
  int i,j;
  
  pi=blockIdx.x*blockDim.x+threadIdx.x;

cuda_cmap[pi]--;

}

__global__ void findc4(int *cuda_match,int *real,int *cmap,int *cnedges,int *xadj)
{

  int  pi;
  int ii,iii;
  int i,j,u;
  
  pi=blockIdx.x*blockDim.x+threadIdx.x;

  if(pi>cuda_match[pi])
  cmap[pi]=cmap[cuda_match[pi]];

  else
  {
    u=cuda_match[pi];
    real[cmap[pi]]=pi;

    if(u!=pi)
    cnedges[cmap[pi]]=(xadj[pi+1]-xadj[pi])+(xadj[u+1]-xadj[u]);
    else
    cnedges[cmap[pi]]=(xadj[pi+1]-xadj[pi]);
  }



}




int Match_SHEM(ctrl_t *ctrl, graph_t *graph,int level)
{
     int nvtxs=graph->nvtxs;


     cudaDeviceSynchronize();
     initcuda_match<<<nvtxs,1>>>(graph->cuda_match);
     cudaDeviceSynchronize();
     

     cudaMemcpy(  graph->cuda_maxvwgt, ctrl->maxvwgt, sizeof(int), cudaMemcpyHostToDevice); 
    // printf("%d\n",level);
     t_cmk1_start=clock(); 
     cudaDeviceSynchronize();
    //  if(graph->nvtxs>1024)
      cuda_shem<<<1024,1>>>(graph->cuda_nvtxs,graph->cuda_match,graph->cuda_xadj,\
      graph-> cuda_vwgt,graph->cuda_adjwgt,graph->cuda_adjncy,graph->cuda_maxvwgt);
    //  else
   //   cuda_shem1<<<1 ,1>>>(graph->cuda_nvtxs,graph->cuda_match,graph->cuda_xadj,\
        graph-> cuda_vwgt,graph->cuda_adjwgt,graph->cuda_adjncy,graph->cuda_maxvwgt);

      cudaDeviceSynchronize();
     t_cmk1_end=clock(); 
     duration_cmk1 += (double)(t_cmk1_end - t_cmk1_start) / CLOCKS_PER_SEC/test_time; 

 t_cmk2_start=clock(); 
      cudaDeviceSynchronize();
      cuda_cleanv<<<nvtxs,1>>>(graph->cuda_match,graph->cuda_s);
      cudaDeviceSynchronize();
 t_cmk2_end=clock(); 
      duration_cmk2 += (double)(t_cmk2_end - t_cmk2_start) / CLOCKS_PER_SEC/test_time; 


       int cnvtxs=0;
       cudaMemcpy(  graph->cuda_cnvtxs, &cnvtxs,  sizeof(int), cudaMemcpyHostToDevice);

       cudaDeviceSynchronize();

 

  int *cuda_temp;
  cudaMalloc((void**)&cuda_temp,  sizeof(int));

  t_cmk3_start=clock(); 
  cudaDeviceSynchronize();
   findc1<<<nvtxs,1>>>(graph->cuda_match,graph->cuda_cmap,graph->cuda_cnvtxs,graph->cuda_s);
   cudaDeviceSynchronize();
   
   t_cmk3_end=clock(); 
   duration_cmk3 += (double)(t_cmk3_end - t_cmk3_start) / CLOCKS_PER_SEC/test_time; 

   t_cmk4_start=clock(); 
   cudaDeviceSynchronize();
   findc2<<<nvtxs,1>>>(graph->cuda_match,graph->cuda_cmap,graph->cuda_cnvtxs,graph->cuda_s);
   cudaDeviceSynchronize();
   t_cmk4_end=clock(); 
   duration_cmk4 += (double)(t_cmk4_end - t_cmk4_start) / CLOCKS_PER_SEC/test_time; 

   thrust::device_ptr<int> cscan = thrust::device_pointer_cast<int>(graph->cuda_cmap);


   t_cmk5_start=clock(); 
   cudaDeviceSynchronize();
   thrust::inclusive_scan(cscan,cscan+nvtxs,cscan);
   cudaDeviceSynchronize();

   t_cmk5_end=clock(); 
   duration_cmk5 += (double)(t_cmk5_end - t_cmk5_start) / CLOCKS_PER_SEC/test_time; 
 

   t_cmk6_start=clock(); 
   cudaDeviceSynchronize();
   findc2_5<<<1,1>>>(cuda_temp,graph->cuda_cmap,graph->cuda_nvtxs);
   cudaDeviceSynchronize();

   t_cmk6_end=clock(); 
   duration_cmk6 += (double)(t_cmk6_end - t_cmk6_start) / CLOCKS_PER_SEC/test_time; 

   cudaMemcpy(  &cnvtxs,cuda_temp,  sizeof(int), cudaMemcpyDeviceToHost);

   cudaFree(cuda_temp);
   int *cpu_scan_edge=(int *)malloc(sizeof(int)*(cnvtxs+1));
      cudaMalloc((void**)&graph->cuda_real_nvtxs,  cnvtxs*sizeof(int));
      cudaMalloc((void**)&graph->cuda_real_edge,  (cnvtxs+1)*sizeof(int));

      t_cmk7_start=clock(); 
      cudaDeviceSynchronize();
   findc3<<<nvtxs,1>>>(graph->cuda_match,graph->cuda_cmap,graph->cuda_cnvtxs,graph->cuda_s);
   cudaDeviceSynchronize();

   t_cmk7_end=clock(); 
   duration_cmk7 += (double)(t_cmk7_end - t_cmk7_start) / CLOCKS_PER_SEC/test_time; 


   t_cmk8_start=clock(); 
   cudaDeviceSynchronize();
   findc4<<<nvtxs,1>>>(graph->cuda_match,graph->cuda_real_nvtxs,graph->cuda_cmap,graph->cuda_real_edge,graph->cuda_xadj);
   cudaDeviceSynchronize();

   t_cmk8_end=clock(); 
   duration_cmk8 += (double)(t_cmk8_end - t_cmk8_start) / CLOCKS_PER_SEC/test_time; 

     cudaMemcpy(  cpu_scan_edge, graph->cuda_real_edge,  (cnvtxs+1)*sizeof(int), cudaMemcpyDeviceToHost);


      CreateCoarseGraph(ctrl, graph, cnvtxs, level,cpu_scan_edge);
      
     return cnvtxs;
     
}


void Memcpycoarsen_initinfo(graph_t *graph)
{
  int nvtxs=graph->nvtxs;
  int nedges=graph->nedges;

  graph->vwgt=(int *)malloc(sizeof(int)*nvtxs); 
  graph->adjncy=(int *)malloc(sizeof(int)*nedges);
  graph->adjwgt=(int *)malloc(sizeof(int)*nedges);

  t_gtoc_start=clock(); 

  cudaDeviceSynchronize();
  cudaMemcpy(  graph->xadj, graph->cuda_xadj , (nvtxs+1)*sizeof(int), cudaMemcpyDeviceToHost);
  cudaMemcpy(  graph->vwgt, graph->cuda_vwgt , (nvtxs)*sizeof(int), cudaMemcpyDeviceToHost);
  cudaMemcpy(  graph->adjncy, graph->cuda_adjncy , (nedges)*sizeof(int), cudaMemcpyDeviceToHost);
  cudaMemcpy(  graph->adjwgt, graph->cuda_adjwgt , (nedges)*sizeof(int), cudaMemcpyDeviceToHost);
  cudaDeviceSynchronize();

  t_gtoc_end=clock(); 
  duration_gtoc += (double)(t_gtoc_end - t_gtoc_start) / CLOCKS_PER_SEC/test_time; 
}


graph_t *CoarsenGraph(ctrl_t *ctrl, graph_t *graph)
{
 
  int i, level=0;

  ctrl->maxvwgt[0] = 1.5*graph->tvwgt[0]/ctrl->CoarsenTo;
/*
  printf("%d\n",ctrl->CoarsenTo);

printf("%d\n",graph->nvtxs);
*/
  do{

          if(level!=0)
          {Malloccoarseninfo(ctrl,graph);}

          Match_SHEM(ctrl,graph,level);
          graph = graph->coarser;
          level++;
       //   printf("%d  %d %d %d %d\n",graph->nvtxs,  ctrl->CoarsenTo ,int(0.75*graph->finer->nvtxs),graph->nedges,graph->nvtxs/2);
          fflush(NULL);
        }while (graph->nvtxs > ctrl->CoarsenTo && \
        graph->nvtxs < 0.75*graph->finer->nvtxs &&\
        graph->nedges > graph->nvtxs/2  );




  //bb_segsort(graph->cuda_adjwgt, graph->cuda_adjncy, graph->nedges, graph->cuda_xadj, graph->nvtxs);


   Memcpycoarsen_initinfo(graph);  



  // printf("%d",graph->nvtxs);

  return graph;
 
}



graph_t *CoarsenGraph1(ctrl_t *ctrl, graph_t *graph)
{
 
  int i, eqewgts, level=1;

  for (eqewgts=1, i=1; i<graph->nedges; i++) {
    if (graph->adjwgt[0] != graph->adjwgt[i]) {
      eqewgts = 0;
      break;
    }
  }

    ctrl->maxvwgt[0] = 1.5*graph->tvwgt[0]/ctrl->CoarsenTo;

  
 do
 {
    if (graph->cmap == NULL)
      graph->cmap = (int*)malloc(sizeof(int)*(graph->nvtxs));

          Match_SHEM0 (ctrl,graph,level);

          graph = graph->coarser;
          
    eqewgts = 0;
    level++;

 }while (graph->nvtxs > ctrl->CoarsenTo && 
           graph->nvtxs < 0.75*graph->finer->nvtxs && 
           graph->nedges > graph->nvtxs/2);

  return graph;
}



void Allocate2WayPartitionMemory(ctrl_t *ctrl, graph_t *graph)
{
  int nvtxs;

  nvtxs = graph->nvtxs;
  graph->pwgts  = (int*)malloc(2*sizeof(int));
  graph->where  = (int*)malloc(nvtxs*sizeof(int));
  graph->bndptr = (int*)malloc(nvtxs*sizeof(int));
  graph->bndind = (int*)malloc(nvtxs*sizeof(int));
  graph->id     = (int*)malloc(nvtxs*sizeof(int));
  graph->ed     = (int*)malloc(nvtxs*sizeof(int));
}

void Compute2WayPartitionParams(ctrl_t *ctrl, graph_t *graph)
{
 int i, j, nvtxs, nbnd, mincut, istart, iend, tid, ted, me;
  int *xadj, *vwgt, *adjncy, *adjwgt, *pwgts;
  int *where, *bndptr, *bndind, *id, *ed;

  nvtxs  = graph->nvtxs;

  xadj   = graph->xadj;
  vwgt   = graph->vwgt;
  adjncy = graph->adjncy;
  adjwgt = graph->adjwgt;

  where  = graph->where;
  id     = graph->id;
  ed     = graph->ed;

  pwgts  = iset(2, 0, graph->pwgts);
  bndptr = iset(nvtxs, -1, graph->bndptr);
  bndind = graph->bndind;


    for (i=0; i<nvtxs; i++) {
      pwgts[where[i]] += vwgt[i];
    }



  for (nbnd=0, mincut=0, i=0; i<nvtxs; i++) {
    istart = xadj[i];
    iend   = xadj[i+1];

    me = where[i];
    tid = ted = 0;
    for (j=istart; j<iend; j++) {
      if (me == where[adjncy[j]])
        tid += adjwgt[j];
      else
        ted += adjwgt[j];
    }
    id[i] = tid;
    ed[i] = ted;
  
    if (ted > 0 || istart == iend) {
      BNDInsert(nbnd, bndind, bndptr, i);
      mincut += ted;
    }
  }

  graph->mincut = mincut/2;
  graph->nbnd   = nbnd;
  
}

float ComputeLoadImbalanceDiff(graph_t *graph, int nparts, float *pijbm,
           float *ubvec)
{
  int  j, *pwgts;
  float max, cur;


  pwgts = graph->pwgts;

  max = -1.0;

    for (j=0; j<nparts; j++) {
      cur = pwgts[j]*pijbm[j] - ubvec[0];
      if (cur > max)
        max = cur;
    }
  

  return max;
}

 void  rpqInit(rpq_t *queue, size_t maxnodes)
{
    int i;
  queue->nnodes = 0;
  queue->maxnodes = maxnodes;

  queue->heap    = (rkv_t*)malloc(sizeof(rkv_t)*maxnodes);
  queue->locator=(ssize_t*)malloc(sizeof(ssize_t)*maxnodes);
  for(i=0;i<maxnodes;i++)
  {
      queue->locator[i]=-1;
  }
}

rpq_t *rpqCreate(size_t maxnodes)
{
  rpq_t *queue; 
  queue = (rpq_t *)malloc(sizeof(rpq_t));
  rpqInit(queue, maxnodes);

  return queue;
}

int    rpqInsert(rpq_t *queue, int node, int key)
{
  ssize_t i, j;
  ssize_t *locator=queue->locator;
  rkv_t *heap=queue->heap;
  i = queue->nnodes++;
  while (i > 0) {
    j = (i-1)>>1;
    if (KEY_LT(key, heap[j].key)) {
      heap[i] = heap[j];
      locator[heap[i].val] = i;
      i = j;
    }
    else
      break;
  }

  heap[i].key   = key;
  heap[i].val   = node;
  locator[node] = i;
  
  return 0;
}

  int     rpqGetTop(rpq_t *queue)
{

  ssize_t i, j;
  ssize_t *locator;
  rkv_t *heap;
  int vtx, node;
  float key;
  if (queue->nnodes == 0)
    return -1;

  queue->nnodes--;

  heap    = queue->heap;
  locator = queue->locator;

  vtx = heap[0].val;
  locator[vtx] = -1;

  if ((i = queue->nnodes) > 0) {
    key  = heap[i].key;
    node = heap[i].val;
    i = 0;
    while ((j=2*i+1) < queue->nnodes) {
      if (KEY_LT(heap[j].key, key)) {
        if (j+1 < queue->nnodes && KEY_LT(heap[j+1].key, heap[j].key))
          j = j+1;
        heap[i] = heap[j];
        locator[heap[i].val] = i;
        i = j;
      }
      else if (j+1 < queue->nnodes && KEY_LT(heap[j+1].key, key)) {
        j = j+1;
        heap[i] = heap[j];
        locator[heap[i].val] = i;
        i = j;
      }
      else
        break;
    }

    heap[i].key   = key;
    heap[i].val   = node;
    locator[node] = i;
  }
  return vtx;
}

  int    rpqDelete(rpq_t *queue, int node)
{
   
  ssize_t i, j, nnodes;
  float newkey, oldkey;
  ssize_t *locator=queue->locator;
  rkv_t *heap=queue->heap;
  i = locator[node];
  locator[node] = -1;

  if (--queue->nnodes > 0 && heap[queue->nnodes].val != node) {
    node   = heap[queue->nnodes].val;
    newkey = heap[queue->nnodes].key;
    oldkey = heap[i].key;

    if (KEY_LT(newkey, oldkey)) { 
      while (i > 0) {
        j = (i-1)>>1;
        if (KEY_LT(newkey, heap[j].key)) {
          heap[i] = heap[j];
          locator[heap[i].val] = i;
          i = j;
        }
        else
          break;
      }
    }
    else { 
      nnodes = queue->nnodes;
      while ((j=(i<<1)+1) < nnodes) {
        if (KEY_LT(heap[j].key, newkey)) {
          if (j+1 < nnodes && KEY_LT(heap[j+1].key, heap[j].key))
            j++;
          heap[i] = heap[j];
          locator[heap[i].val] = i;
          i = j;
        }
        else if (j+1 < nnodes && KEY_LT(heap[j+1].key, newkey)) {
          j++;
          heap[i] = heap[j];
          locator[heap[i].val] = i;
          i = j;
        }
        else
          break;
      }
    }

    heap[i].key   = newkey;
    heap[i].val   = node;
    locator[node] = i;
  }
  return 0;
}
  void   rpqUpdate(rpq_t *queue, int node, int newkey)
{
 
  ssize_t i, j, nnodes;
  float oldkey;
  ssize_t *locator=queue->locator;
  rkv_t *heap=queue->heap;

  oldkey = heap[locator[node]].key;
  i = locator[node];

  if (KEY_LT(newkey, oldkey)) { 
    while (i > 0) {
      j = (i-1)>>1;
      if (KEY_LT(newkey, heap[j].key)) {
        heap[i] = heap[j];
        locator[heap[i].val] = i;
        i = j;
      }
      else
        break;
    }
  }
  else { 
    nnodes = queue->nnodes;
    while ((j=(i<<1)+1) < nnodes) {
      if (KEY_LT(heap[j].key, newkey)) {
        if (j+1 < nnodes && KEY_LT(heap[j+1].key, heap[j].key))
          j++;
        heap[i] = heap[j];
        locator[heap[i].val] = i;
        i = j;
      }
      else if (j+1 < nnodes && KEY_LT(heap[j+1].key, newkey)) {
        j++;
        heap[i] = heap[j];
        locator[heap[i].val] = i;
        i = j;
      }
      else
        break;
    }
  }

  heap[i].key   = newkey;
  heap[i].val   = node;
  locator[node] = i;
  return;
}

void gk_gkmcoreDel(gk_mcore_t *mcore, void *ptr)
{
  int i;

  for (i=mcore->cmop-1; i>=0; i--) {
    if (mcore->mops[i].type == 1)
      exit(0);

    if (mcore->mops[i].ptr == ptr) {
      if (mcore->mops[i].type != 3)
        exit(0);

      mcore->cur_hallocs -= mcore->mops[i].nbytes;
      mcore->mops[i] = mcore->mops[--mcore->cmop];
      return;
    }
  }

  exit(0);
}

void gk_free(void **ptr1,...)
{
  va_list plist;
  void **ptr;

  if (*ptr1 != NULL) {
    free(*ptr1);

    if (gkmcore != NULL) gk_gkmcoreDel(gkmcore, *ptr1);
  }
  *ptr1 = NULL;

  va_start(plist, ptr1);
  while ((ptr = va_arg(plist, void **)) != LTERM) {
    if (*ptr != NULL) {
      free(*ptr);


      if (gkmcore != NULL) gk_gkmcoreDel(gkmcore, *ptr);
    }
    *ptr = NULL;
  }
  va_end(plist);
}       

  void   rpqFree(rpq_t *queue)
{
  if (queue == NULL) return;
  gk_free((void **)&queue->heap, &queue->locator, LTERM);
  queue->maxnodes = 0;
}

 void   rpqDestroy(rpq_t *queue)
  {
  if (queue == NULL) return;
  rpqFree(queue);
  gk_free((void **)&queue, LTERM);
}

 void   rpqReset(rpq_t *queue)
  {
  ssize_t i;
  ssize_t *locator=queue->locator;
  rkv_t *heap=queue->heap;

  for (i=queue->nnodes-1; i>=0; i--)
    locator[heap[i].val] = -1;
  queue->nnodes = 0;
}

void Bnd2WayBalance(ctrl_t *ctrl, graph_t *graph, float *ntpwgts)
{
  int i, ii, j, k, kwgt, nvtxs, nbnd, nswaps, from, to, pass, me, tmp;
  int *xadj, *vwgt, *adjncy, *adjwgt, *where, *id, *ed, *bndptr, *bndind, *pwgts;
  int *moved, *perm;
  rpq_t *queue;
  int higain, mincut, mindiff;
  int tpwgts[2];



  nvtxs  = graph->nvtxs;
  xadj   = graph->xadj;
  vwgt   = graph->vwgt;
  adjncy = graph->adjncy;
  adjwgt = graph->adjwgt;
  where  = graph->where;
  id     = graph->id;
  ed     = graph->ed;
  pwgts  = graph->pwgts;
  bndptr = graph->bndptr;
  bndind = graph->bndind;

  moved = iwspacemalloc(ctrl, nvtxs);
  perm  = iwspacemalloc(ctrl, nvtxs);

  tpwgts[0] = graph->tvwgt[0]*ntpwgts[0];
  tpwgts[1] = graph->tvwgt[0] - tpwgts[0];
  mindiff   = abs(tpwgts[0]-pwgts[0]);
  from      = (pwgts[0] < tpwgts[0] ? 1 : 0);
  to        = (from+1)%2;

  queue = rpqCreate(nvtxs);

  iset(nvtxs, -1, moved);
  nbnd = graph->nbnd;
  irandArrayPermute(nbnd, perm, nbnd/5, 1);
  for (ii=0; ii<nbnd; ii++) {
    i = perm[ii];
    if (where[bndind[i]] == from && vwgt[bndind[i]] <= mindiff)
      rpqInsert(queue, bndind[i], ed[bndind[i]]-id[bndind[i]]);
  }

  mincut = graph->mincut;
  
  for (nswaps=0; nswaps<nvtxs; nswaps++) 
  {
    if ((higain = rpqGetTop(queue)) == -1)
      break;
    if (pwgts[to]+vwgt[higain] > tpwgts[to])
      break;
    mincut -= (ed[higain]-id[higain]);
    INC_DEC(pwgts[to], pwgts[from], vwgt[higain]);
    where[higain] = to;
    moved[higain] = nswaps;
    gk_SWAP(id[higain], ed[higain], tmp);
    if (ed[higain] == 0 && xadj[higain] < xadj[higain+1]) 
      BNDDelete(nbnd, bndind,  bndptr, higain);
    for (j=xadj[higain]; j<xadj[higain+1]; j++) {
      k = adjncy[j];
      kwgt = (to == where[k] ? adjwgt[j] : -adjwgt[j]);
      INC_DEC(id[k], ed[k], kwgt);


      if (bndptr[k] != -1) { 
        if (ed[k] == 0) { 
          BNDDelete(nbnd, bndind, bndptr, k);
          if (moved[k] == -1 && where[k] == from && vwgt[k] <= mindiff) 
            rpqDelete(queue, k);
        }
        else { 
          if (moved[k] == -1 && where[k] == from && vwgt[k] <= mindiff)
            rpqUpdate(queue, k, ed[k]-id[k]);
        }
      }
      else {
        if (ed[k] > 0) {  
          BNDInsert(nbnd, bndind, bndptr, k);
          if (moved[k] == -1 && where[k] == from && vwgt[k] <= mindiff) 
            rpqInsert(queue, k, ed[k]-id[k]);
        }
      }
    }
  }


  graph->mincut = mincut;
  graph->nbnd   = nbnd;

  rpqDestroy(queue);

}

void Balance2Way(ctrl_t *ctrl, graph_t *graph, float *ntpwgts)
{
  if (ComputeLoadImbalanceDiff(graph, 2, ctrl->pijbm, ctrl->ubfactors) <= 0) 
    return;

    if (abs(ntpwgts[0]*graph->tvwgt[0]-graph->pwgts[0]) < 3*graph->tvwgt[0]/graph->nvtxs)
      return;

      Bnd2WayBalance(ctrl, graph, ntpwgts);

  

}


void FM_2WayCutRefine(ctrl_t *ctrl, graph_t *graph, float *ntpwgts, int niter)
{

  int i, ii, j, k, kwgt, nvtxs, nbnd, nswaps, from, to, pass, me, limit, tmp;
  int *xadj, *vwgt, *adjncy, *adjwgt, *where, *id, *ed, *bndptr, *bndind, *pwgts;
  int *moved, *swaps, *perm;
  rpq_t *queues[2];
  int higain, mincut, mindiff, origdiff, initcut, newcut, mincutorder, avgvwgt;
  int tpwgts[2];

  nvtxs  = graph->nvtxs;
  xadj   = graph->xadj;
  vwgt   = graph->vwgt;
  adjncy = graph->adjncy;
  adjwgt = graph->adjwgt;
  where  = graph->where;
  id     = graph->id;
  ed     = graph->ed;
  pwgts  = graph->pwgts;
  bndptr = graph->bndptr;
  bndind = graph->bndind;

  moved = iwspacemalloc(ctrl, nvtxs);
  swaps = iwspacemalloc(ctrl, nvtxs);
  perm  = iwspacemalloc(ctrl, nvtxs);
  tpwgts[0] = graph->tvwgt[0]*ntpwgts[0];
  tpwgts[1] = graph->tvwgt[0]-tpwgts[0];
  limit   = gk_min(gk_max(0.01*nvtxs, 15), 100);
  avgvwgt = gk_min((pwgts[0]+pwgts[1])/20, 2*(pwgts[0]+pwgts[1])/nvtxs);
  queues[0] = rpqCreate(nvtxs);
  queues[1] = rpqCreate(nvtxs);
  origdiff = abs(tpwgts[0]-pwgts[0]);
  iset(nvtxs, -1, moved);
  for (pass=0; pass<niter; pass++) 
  { 
  
    rpqReset(queues[0]);
    rpqReset(queues[1]);
 
    mincutorder = -1;
    newcut = mincut = initcut = graph->mincut;
    mindiff = abs(tpwgts[0]-pwgts[0]);

    nbnd = graph->nbnd;
    irandArrayPermute(nbnd, perm, nbnd, 1);  

    for (ii=0; ii<nbnd; ii++) {
      i = perm[ii];
      rpqInsert(queues[where[bndind[i]]], bndind[i], ed[bndind[i]]-id[bndind[i]]);
    }   
    
    for (nswaps=0; nswaps<nvtxs; nswaps++) {
      from = (tpwgts[0]-pwgts[0] < tpwgts[1]-pwgts[1] ? 0 : 1);
      to = (from+1)%2;

      if ((higain = rpqGetTop(queues[from])) == -1)
        break;

      newcut -= (ed[higain]-id[higain]);
      INC_DEC(pwgts[to], pwgts[from], vwgt[higain]);

      if ((newcut < mincut && abs(tpwgts[0]-pwgts[0]) <= origdiff+avgvwgt) || 
          (newcut == mincut && abs(tpwgts[0]-pwgts[0]) < mindiff)) {
        mincut  = newcut;
        mindiff = abs(tpwgts[0]-pwgts[0]);
        mincutorder = nswaps;
      }
      else if (nswaps-mincutorder > limit) { 
        newcut += (ed[higain]-id[higain]);
        INC_DEC(pwgts[from], pwgts[to], vwgt[higain]);
        break;
      }

      where[higain] = to;
      moved[higain] = nswaps;
      swaps[nswaps] = higain;
      gk_SWAP(id[higain], ed[higain], tmp);
      if (ed[higain] == 0 && xadj[higain] < xadj[higain+1]) 
        BNDDelete(nbnd, bndind,  bndptr, higain);

      for (j=xadj[higain]; j<xadj[higain+1]; j++) {
        k = adjncy[j];
        kwgt = (to == where[k] ? adjwgt[j] : -adjwgt[j]);
        INC_DEC(id[k], ed[k], kwgt);
        if (bndptr[k] != -1) { 
          if (ed[k] == 0) { 
            BNDDelete(nbnd, bndind, bndptr, k);
            if (moved[k] == -1)  
              rpqDelete(queues[where[k]], k);
          }
          else { 
            if (moved[k] == -1) 
              rpqUpdate(queues[where[k]], k, ed[k]-id[k]);
          }
        }
        else {
          if (ed[k] > 0) {  
            BNDInsert(nbnd, bndind, bndptr, k);
            if (moved[k] == -1) 
              rpqInsert(queues[where[k]], k, ed[k]-id[k]);
          }
        }
      }
    }
    for (i=0; i<nswaps; i++)
      moved[swaps[i]] = -1;  
    for (nswaps--; nswaps>mincutorder; nswaps--) {
      higain = swaps[nswaps];
      to = where[higain] = (where[higain]+1)%2;
      gk_SWAP(id[higain], ed[higain], tmp);
      if (ed[higain] == 0 && bndptr[higain] != -1 && xadj[higain] < xadj[higain+1])
        BNDDelete(nbnd, bndind,  bndptr, higain);
      else if (ed[higain] > 0 && bndptr[higain] == -1)
        BNDInsert(nbnd, bndind,  bndptr, higain);
      INC_DEC(pwgts[to], pwgts[(to+1)%2], vwgt[higain]);
      for (j=xadj[higain]; j<xadj[higain+1]; j++) {
        k = adjncy[j];

        kwgt = (to == where[k] ? adjwgt[j] : -adjwgt[j]);
        INC_DEC(id[k], ed[k], kwgt);

        if (bndptr[k] != -1 && ed[k] == 0)
          BNDDelete(nbnd, bndind, bndptr, k);
        if (bndptr[k] == -1 && ed[k] > 0)
          BNDInsert(nbnd, bndind, bndptr, k);
      }
    }

    graph->mincut = mincut;
    graph->nbnd   = nbnd;

    if (mincutorder <= 0 || mincut == initcut)
      break;
  }
   

  rpqDestroy(queues[0]);
  rpqDestroy(queues[1]);

}



void FM_2WayRefine(ctrl_t *ctrl, graph_t *graph, float *ntpwgts, float niter)
{
    FM_2WayCutRefine(ctrl, graph, ntpwgts, niter);
}


void GrowBisection(ctrl_t *ctrl, graph_t *graph, float *ntpwgts, 
         int niparts)
{
  int i, j, k, nvtxs, drain, nleft, first, last, 
        pwgts[2], oneminpwgt, onemaxpwgt, 
        from, me, bestcut=0, icut, mincut, inbfs;
  int *xadj, *vwgt, *adjncy, *adjwgt, *where;
  int *queue, *touched, *gain, *bestwhere;
  

  nvtxs  = graph->nvtxs;
  xadj   = graph->xadj;
  vwgt   = graph->vwgt;
  adjncy = graph->adjncy;
  adjwgt = graph->adjwgt;
  Allocate2WayPartitionMemory(ctrl, graph);
  where = graph->where;
  bestwhere = iwspacemalloc(ctrl, nvtxs);
  queue     = iwspacemalloc(ctrl, nvtxs);
  touched   = iwspacemalloc(ctrl, nvtxs);
  onemaxpwgt = ctrl->ubfactors[0]*graph->tvwgt[0]*ntpwgts[1];
  oneminpwgt = (1.0/ctrl->ubfactors[0])*graph->tvwgt[0]*ntpwgts[1]; 
  
  
  for (inbfs=0; inbfs<niparts; inbfs++) 
  {
    iset(nvtxs, 1, where);

    iset(nvtxs, 0, touched);

    pwgts[1] = graph->tvwgt[0];
    pwgts[0] = 0;
    queue[0] = irandInRange(nvtxs);
    touched[queue[0]] = 1;
    first = 0; 
    last  = 1;
    nleft = nvtxs-1;
    drain = 0;
    for (;;)        
    {
      if (first == last) 
      { 

        if (nleft == 0 || drain)
          break;

        k = irandInRange(nleft);

        for (i=0; i<nvtxs; i++) 
        {
          if (touched[i] == 0) {
            if (k == 0)
              break;
            else
              k--;
          }
        }

        queue[0]   = i;
        touched[i] = 1;
        first      = 0; 
        last       = 1;
        nleft--;

      }

      i = queue[first++];
      if (pwgts[0] > 0 && pwgts[1]-vwgt[i] < oneminpwgt) {
        drain = 1;
        continue;
      }

      where[i] = 0;
      INC_DEC(pwgts[0], pwgts[1], vwgt[i]);
      if (pwgts[1] <= onemaxpwgt)
        break;

      drain = 0;
   
      for (j=xadj[i]; j<xadj[i+1]; j++) {
        k = adjncy[j];
        if (touched[k] == 0) {
          queue[last++] = k;
          touched[k] = 1;
          nleft--;
        }
      }
    }
/*
    if (pwgts[1] == 0) 
      where[irandInRange(nvtxs)] = 1;
    if (pwgts[0] == 0) 
      where[irandInRange(nvtxs)] = 0; 
*/

    Compute2WayPartitionParams(ctrl, graph);
    Balance2Way(ctrl, graph, ntpwgts);
    FM_2WayRefine(ctrl, graph, ntpwgts, ctrl->niter);

    if (inbfs == 0 || bestcut > graph->mincut) {
      bestcut = graph->mincut;
      icopy(nvtxs, where, bestwhere);
      if (bestcut == 0)
        break;
    }
  }

  graph->mincut = bestcut;
  icopy(nvtxs, bestwhere, where);


}

void Init2WayPartition(ctrl_t *ctrl, graph_t *graph, float *ntpwgts, 
         int niparts) 
{
    GrowBisection(ctrl, graph, ntpwgts, niparts);
}

void FreeRData(graph_t *graph) 
{

  if ((void *)graph->ckrinfo == (void *)graph->vkrinfo)
    graph->ckrinfo = NULL;



  gk_free((void **)&graph->where, &graph->pwgts, &graph->id, &graph->ed, 
      &graph->bndptr, &graph->bndind,  &graph->ckrinfo, 
      &graph->vkrinfo, LTERM);
}

void FreeGraph(graph_t **r_graph) 
{
  graph_t *graph;
  graph = *r_graph;

  if (graph->free_xadj)
    gk_free((void **)&graph->xadj, LTERM);
  if (graph->free_vwgt)
    gk_free((void **)&graph->vwgt, LTERM);
  if (graph->free_vsize)
    gk_free((void **)&graph->vsize, LTERM);
  if (graph->free_adjncy)
    gk_free((void **)&graph->adjncy, LTERM);
  if (graph->free_adjwgt)
    gk_free((void **)&graph->adjwgt, LTERM);

  FreeRData(graph);

  gk_free((void **)&graph->tvwgt, &graph->invtvwgt, &graph->label, 
      &graph->cmap, &graph, LTERM);

  *r_graph = NULL;
}

void Project2WayPartition(ctrl_t *ctrl, graph_t *graph)
{

  int i, j, istart, iend, nvtxs, nbnd, me, tid, ted;
  int *xadj, *adjncy, *adjwgt;
  int *cmap, *where, *bndptr, *bndind;
  int *cwhere, *cbndptr;
  int *id, *ed;
  graph_t *cgraph;

  Allocate2WayPartitionMemory(ctrl, graph);

  cgraph  = graph->coarser;

  cwhere  = cgraph->where;
  
  cbndptr = cgraph->bndptr;

  nvtxs   = graph->nvtxs;
  cmap    = graph->cmap;
  xadj    = graph->xadj;
  adjncy  = graph->adjncy;
  adjwgt  = graph->adjwgt;

  where  = graph->where;
  id     = graph->id;
  ed     = graph->ed;

  bndptr = iset(nvtxs, -1, graph->bndptr);
  bndind = graph->bndind;

  for (i=0; i<nvtxs; i++) {
    j = cmap[i];
    where[i] = cwhere[j];
    cmap[i]  = cbndptr[j];
  }


  
  for (nbnd=0, i=0; i<nvtxs; i++) {
    istart = xadj[i];
    iend   = xadj[i+1];
  
    tid = ted = 0;
    if (cmap[i] == -1) { 
    
      for (j=istart; j<iend; j++)
        tid += adjwgt[j];
    }
    else { 
      me = where[i];
      for (j=istart; j<iend; j++) {
        if (me == where[adjncy[j]])
          tid += adjwgt[j];
        else
          ted += adjwgt[j];
      }
    }
    id[i] = tid;
    ed[i] = ted;

    if (ted > 0 || istart == iend) 
      BNDInsert(nbnd, bndind, bndptr, i);
  }
  graph->mincut = cgraph->mincut;
  graph->nbnd   = nbnd;

  icopy(2, cgraph->pwgts, graph->pwgts);

  FreeGraph(&graph->coarser);
  graph->coarser = NULL;
  
}

void Refine2Way(ctrl_t *ctrl, graph_t *orggraph, graph_t *graph, float *tpwgts)
{
 

  Compute2WayPartitionParams(ctrl, graph);

  for (;;) {
    Balance2Way(ctrl, graph, tpwgts);

    FM_2WayRefine(ctrl, graph, tpwgts, ctrl->niter); 
  
    if (graph == orggraph)
      break;
  
    graph = graph->finer;
    Project2WayPartition(ctrl, graph);


  }

   
}

int MultilevelBisect(ctrl_t *ctrl, graph_t *graph, float *tpwgts)
{

  int i, niparts, bestobj=0, curobj=0, *bestwhere=NULL;
  graph_t *cgraph;
  float bestbal=0.0, curbal=0.0;
  Setup2WayBalMultipliers(ctrl, graph, tpwgts);
   cgraph = CoarsenGraph1(ctrl, graph);

   niparts = 10;
   Init2WayPartition(ctrl, cgraph, tpwgts, niparts);
    
   Refine2Way(ctrl, graph, cgraph, tpwgts);

   
    curobj = graph->mincut;
    curbal = ComputeLoadImbalanceDiff(graph, 2, ctrl->pijbm, ctrl->ubfactors);

      bestobj = curobj;
      bestbal = curbal;

  if (bestobj != curobj) {
    icopy(graph->nvtxs, bestwhere, graph->where);
    Compute2WayPartitionParams(ctrl, graph);
  
  }
  


  return bestobj;
}



graph_t *SetupSplitGraph(graph_t *graph, int snvtxs, int snedges)
{
  graph_t *sgraph;

  sgraph = CreateGraph();
  sgraph->nvtxs  = snvtxs;
  sgraph->nedges = snedges;

  sgraph->xadj        = (int*)malloc(sizeof(int)*(snvtxs+1));
  sgraph->vwgt        = (int*)malloc(sizeof(int)*(snvtxs+1));
  sgraph->adjncy      = (int*)malloc(sizeof(int)*(snedges));
  sgraph->adjwgt      = (int*)malloc(sizeof(int)*(snedges));
  sgraph->label	      = (int*)malloc(sizeof(int)*(snvtxs));
  sgraph->tvwgt       = (int*)malloc(sizeof(int));
  sgraph->invtvwgt    = (float*)malloc(sizeof(float));

  if (graph->vsize)
    sgraph->vsize     = (int*)malloc(sizeof(int)*(snvtxs));

  return sgraph;
}


void SplitGraphPart(ctrl_t *ctrl, graph_t *graph, graph_t **r_lgraph, 
         graph_t **r_rgraph)
{

  int i, j, k, l, istart, iend, mypart, nvtxs, ncon, snvtxs[2], snedges[2];
  int *xadj, *vwgt, *adjncy, *adjwgt, *label, *where, *bndptr;
  int *sxadj[2], *svwgt[2], *sadjncy[2], *sadjwgt[2], *slabel[2];
  int *rename;
  int *auxadjncy, *auxadjwgt;
  graph_t *lgraph, *rgraph;


  nvtxs   = graph->nvtxs;
  xadj    = graph->xadj;
  vwgt    = graph->vwgt;
  adjncy  = graph->adjncy;
  adjwgt  = graph->adjwgt;
  label   = graph->label;
  where   = graph->where;
  bndptr  = graph->bndptr;


  rename = iwspacemalloc(ctrl, nvtxs);
  
  snvtxs[0] = snvtxs[1] = snedges[0] = snedges[1] = 0;
  
  for (i=0; i<nvtxs; i++) {
    k = where[i];
    rename[i] = snvtxs[k]++;
    snedges[k] += xadj[i+1]-xadj[i];
  }

  lgraph      = SetupSplitGraph(graph, snvtxs[0], snedges[0]);
  sxadj[0]    = lgraph->xadj;
  svwgt[0]    = lgraph->vwgt;
  sadjncy[0]  = lgraph->adjncy; 	
  sadjwgt[0]  = lgraph->adjwgt; 
  slabel[0]   = lgraph->label;

  rgraph      = SetupSplitGraph(graph, snvtxs[1], snedges[1]);
  sxadj[1]    = rgraph->xadj;
  svwgt[1]    = rgraph->vwgt;
  sadjncy[1]  = rgraph->adjncy; 	
  sadjwgt[1]  = rgraph->adjwgt; 
  slabel[1]   = rgraph->label;

  snvtxs[0] = snvtxs[1] = snedges[0] = snedges[1] = 0;
  sxadj[0][0] = sxadj[1][0] = 0;

  for (i=0; i<nvtxs; i++) 
  {
    mypart = where[i];
    istart = xadj[i];
    iend = xadj[i+1];
    if (bndptr[i] == -1) 
    { 
      auxadjncy = sadjncy[mypart] + snedges[mypart] - istart;
      auxadjwgt = sadjwgt[mypart] + snedges[mypart] - istart;

      for(j=istart; j<iend; j++) 
      {
        auxadjncy[j] = adjncy[j];
        auxadjwgt[j] = adjwgt[j]; 
      }
      snedges[mypart] += iend-istart;
    }
    else 
    {
      auxadjncy = sadjncy[mypart];
      auxadjwgt = sadjwgt[mypart];
      l = snedges[mypart];
      for (j=istart; j<iend; j++) 
      {
        k = adjncy[j];
        if (where[k] == mypart) {
          auxadjncy[l] = k;
          auxadjwgt[l++] = adjwgt[j]; 
        }
      }
      snedges[mypart] = l;
    }


    
      svwgt[mypart][snvtxs[mypart]] = vwgt[i];

    slabel[mypart][snvtxs[mypart]]   = label[i];
    sxadj[mypart][++snvtxs[mypart]]  = snedges[mypart];
  }

  for (mypart=0; mypart<2; mypart++) {
    iend = sxadj[mypart][snvtxs[mypart]];
    auxadjncy = sadjncy[mypart];
    for (i=0; i<iend; i++) 
      auxadjncy[i] = rename[auxadjncy[i]];
  }

  lgraph->nedges = snedges[0];
  rgraph->nedges = snedges[1];

  SetupGraph_tvwgt(lgraph);
  SetupGraph_tvwgt(rgraph);

  *r_lgraph = lgraph;
  *r_rgraph = rgraph;

}

int MlevelRecursiveBisection(ctrl_t *ctrl, graph_t *graph, int nparts, 
          int *part, float *tpwgts, int fpart)
{
  int i, j, nvtxs, ncon, objval;
  int *label, *where;
  graph_t *lgraph, *rgraph;
  float wsum, *tpwgts2;
 // printf("nvtxs=%d \n",graph->nvtxs);
  if(graph->nvtxs==0)
  {
    printf("****You are trying to partition too many parts!****\n");
    fflush(NULL);
    return 0;
  }
  nvtxs = graph->nvtxs;
  tpwgts2 = rwspacemalloc(ctrl);
    tpwgts2[0]      = rsum((nparts>>1), tpwgts);
   tpwgts2[1] = 1.0 - tpwgts2[0];
    objval = MultilevelBisect(ctrl, graph, tpwgts2);
    
   label = graph->label;
   where = graph->where;

  for (i=0; i<nvtxs; i++)
    part[label[i]] = where[i] + fpart;

  for (i=0; i<nvtxs; i++)
    part[label[i]] = where[i] + fpart;
  if (nparts > 2) 
    SplitGraphPart(ctrl, graph, &lgraph, &rgraph);

  FreeGraph(&graph);

    wsum = rsum((nparts>>1), tpwgts);
     rscale((nparts>>1), 1.0/wsum, tpwgts);
    rscale(nparts-(nparts>>1), 1.0/(1.0-wsum), tpwgts+(nparts>>1));

  if (nparts > 3) {
    objval += MlevelRecursiveBisection(ctrl, lgraph, (nparts>>1), part, tpwgts, fpart);
    objval += MlevelRecursiveBisection(ctrl, rgraph, nparts-(nparts>>1), part, tpwgts+(nparts>>1), fpart+(nparts>>1));
  }
  else if (nparts == 3) {
    FreeGraph(&lgraph);
    objval += MlevelRecursiveBisection(ctrl, rgraph, nparts-(nparts>>1), part, tpwgts+(nparts>>1), fpart+(nparts>>1));
  }
  
    return objval;
}

void SetupKWayBalMultipliers(ctrl_t *ctrl, graph_t *graph)
{
  int i, j;

  for (i=0; i<ctrl->nparts; i++) {
    for (j=0; j<1; j++)
      ctrl->pijbm[i+j] = graph->invtvwgt[j]/ctrl->tpwgts[i+j];
  }
}

int gk_log2(int a)
{
  int i;

  for (i=1; a > 1; i++, a = a>>1);
  return i-1;
}





int PartGraph1(int *nvtxs,  int *xadj, int *adjncy, int *vwgt,int *adjwgt, \
int *nparts, float *tpwgts, float *ubvec, int *objval, int *part)
{
                     graph_t *graph;
                    ctrl_t *ctrl;
                    ctrl = SetupCtrl( *nparts, tpwgts, ubvec);
                    graph = SetupGraph(ctrl, *nvtxs, xadj, adjncy, vwgt, adjwgt);
                    AllocateWorkSpace(ctrl, graph);
                   
                   *objval = MlevelRecursiveBisection(ctrl, graph, *nparts, part, ctrl->tpwgts, 0);
                    return 1;
 }


__global__ void  projectback(int *where,int *cwhere,int *cmap)
{
  int ii,iii,pi;
    
  pi=blockIdx.x*blockDim.x+threadIdx.x;

  where[pi]=cwhere[cmap[pi]];
}



void ProjectKWayPartition(ctrl_t *ctrl, graph_t *graph)
{
            
            int nvtxs=graph->nvtxs;
            graph_t *cgraph; 
            cgraph = graph->coarser;
            t_rpk1_start=clock(); 
            cudaDeviceSynchronize();
            projectback<<<nvtxs,1>>>(graph->cuda_where,cgraph->cuda_where,graph->cuda_cmap);
            cudaDeviceSynchronize();
            t_rpk1_end=clock();
            duration_rpk1 += (double)(t_rpk1_end - t_rpk1_start) / CLOCKS_PER_SEC/test_time; 
}


void InitKWayPartitioning(ctrl_t *ctrl, graph_t *graph)
{
  int i, ntrials, curobj=0, bestobj=0;
  int *bestwhere=NULL;
  float *ubvec=NULL;
  int status;

  graph->where=(int *)malloc(sizeof(int)*graph->nvtxs);
  ubvec = (float*)malloc(sizeof(float));
  ubvec[0] = (float)pow(ctrl->ubfactors[0], 1.0/log(ctrl->nparts));
//printf("%d",graph->nvtxs);
 // printf("ub=%f\n",ubvec[0]);

  PartGraph1(&graph->nvtxs, graph->xadj, graph->adjncy, graph->vwgt, graph->adjwgt, &ctrl->nparts, ctrl->tpwgts, ubvec, \
                           &curobj, graph->where);
  free(ubvec);
  free(bestwhere);

  
}

int IsBalanced(ctrl_t *ctrl, graph_t *graph, float ffactor)
{
  return 
    (ComputeLoadImbalanceDiff(graph, ctrl->nparts, ctrl->pijbm, ctrl->ubfactors) 
         <= ffactor);
}


void cnbrpoolReset(ctrl_t *ctrl)
{
  ctrl->nbrpoolcpos = 0;
}

int cnbrpoolGetNext(ctrl_t *ctrl, int nnbrs)
{
  ctrl->nbrpoolcpos += nnbrs;

  if (ctrl->nbrpoolcpos > ctrl->nbrpoolsize) {
    ctrl->nbrpoolsize += gk_max(10*nnbrs, ctrl->nbrpoolsize/2);

    ctrl->cnbrpool = (cnbr_t *)realloc(ctrl->cnbrpool,  
                          ctrl->nbrpoolsize*sizeof(cnbr_t));
    ctrl->nbrpoolreallocs++;
  }

  return ctrl->nbrpoolcpos - nnbrs;
}




#define UpdateMovedVertexInfoAndBND(i, from, k, to, myrinfo, mynbrs, where, \
            nbnd, bndptr, bndind, bndtype) \
   do { \
     where[i] = to; \
     myrinfo->ed += myrinfo->id-mynbrs[k].ed; \
     gk_SWAP(myrinfo->id, mynbrs[k].ed, j); \
     if (mynbrs[k].ed == 0) \
       mynbrs[k] = mynbrs[--myrinfo->nnbrs]; \
     else \
       mynbrs[k].pid = from; \
     \
     if (bndtype == 1) { \
       if (bndptr[i] != -1 && myrinfo->ed - myrinfo->id < 0) \
         BNDDelete(nbnd, bndind, bndptr, i); \
       if (bndptr[i] == -1 && myrinfo->ed - myrinfo->id >= 0) \
         BNDInsert(nbnd, bndind, bndptr, i); \
     } \
     else { \
       if (bndptr[i] != -1 && myrinfo->ed <= 0) \
         BNDDelete(nbnd, bndind, bndptr, i); \
       if (bndptr[i] == -1 && myrinfo->ed > 0) \
         BNDInsert(nbnd, bndind, bndptr, i); \
     } \
   } while(0) 


#define UpdateAdjacentVertexInfoAndBND(ctrl, vid, adjlen, me, from, to, \
            myrinfo, ewgt, nbnd, bndptr, bndind, bndtype) \
   do { \
     int k; \
     cnbr_t *mynbrs; \
     \
     if (myrinfo->inbr == -1) { \
       myrinfo->inbr  = cnbrpoolGetNext(ctrl, adjlen+1); \
       myrinfo->nnbrs = 0; \
     } \
     \
     mynbrs = ctrl->cnbrpool + myrinfo->inbr; \
     \
     if (me == from) { \
       INC_DEC(myrinfo->ed, myrinfo->id, (ewgt)); \
       if (bndtype == 1) { \
         if (myrinfo->ed-myrinfo->id >= 0 && bndptr[(vid)] == -1) \
           BNDInsert(nbnd, bndind, bndptr, (vid)); \
       } \
       else { \
         if (myrinfo->ed > 0 && bndptr[(vid)] == -1) \
           BNDInsert(nbnd, bndind, bndptr, (vid)); \
       } \
     } \
     else if (me == to) { \
       INC_DEC(myrinfo->id, myrinfo->ed, (ewgt)); \
       if (bndtype == 1) { \
         if (myrinfo->ed-myrinfo->id < 0 && bndptr[(vid)] != -1) \
           BNDDelete(nbnd, bndind, bndptr, (vid)); \
       } \
       else { \
         if (myrinfo->ed <= 0 && bndptr[(vid)] != -1) \
           BNDDelete(nbnd, bndind, bndptr, (vid)); \
       } \
     } \
     \
     if (me != from) { \
       for (k=0; k<myrinfo->nnbrs; k++) { \
         if (mynbrs[k].pid == from) { \
           if (mynbrs[k].ed == (ewgt)) \
             mynbrs[k] = mynbrs[--myrinfo->nnbrs]; \
           else \
             mynbrs[k].ed -= (ewgt); \
           break; \
         } \
       } \
     } \
     \
     if (me != to) { \
       for (k=0; k<myrinfo->nnbrs; k++) { \
         if (mynbrs[k].pid == to) { \
           mynbrs[k].ed += (ewgt); \
           break; \
         } \
       } \
       if (k == myrinfo->nnbrs) { \
         mynbrs[k].pid  = to; \
         mynbrs[k].ed   = (ewgt); \
         myrinfo->nnbrs++; \
       } \
     } \
     \
   } while(0) 



#define UpdateQueueInfo(queue, vstatus, vid, me, from, to, myrinfo, oldnnbrs, \
            nupd, updptr, updind, bndtype) \
   do { \
     float rgain; \
     \
     if (me == to || me == from || oldnnbrs != myrinfo->nnbrs) {  \
       rgain = (myrinfo->nnbrs > 0 ?  \
                1.0*myrinfo->ed/sqrt(myrinfo->nnbrs) : 0.0) - myrinfo->id; \
   \
       if (bndtype == 1) { \
         if (vstatus[(vid)] == 1) { \
           if (myrinfo->ed-myrinfo->id >= 0) \
             rpqUpdate(queue, (vid), rgain); \
           else { \
             rpqDelete(queue, (vid)); \
             vstatus[(vid)] = 3; \
             ListDelete(nupd, updind, updptr, (vid)); \
           } \
         } \
         else if (vstatus[(vid)] == 3 && myrinfo->ed-myrinfo->id >= 0) { \
           rpqInsert(queue, (vid), rgain); \
           vstatus[(vid)] = 1; \
           ListInsert(nupd, updind, updptr, (vid)); \
         } \
       } \
       else { \
         if (vstatus[(vid)] == 1) { \
           if (myrinfo->ed > 0) \
             rpqUpdate(queue, (vid), rgain); \
           else { \
             rpqDelete(queue, (vid)); \
             vstatus[(vid)] = 3; \
             ListDelete(nupd, updind, updptr, (vid)); \
           } \
         } \
         else if (vstatus[(vid)] == 3 && myrinfo->ed > 0) { \
           rpqInsert(queue, (vid), rgain); \
           vstatus[(vid)] = 1; \
           ListInsert(nupd, updind, updptr, (vid)); \
         } \
       } \
     } \
} while(0) 





__global__ void find_real_bnd_info(int *cuda_ppnvtxs,int *cuda_real_bnd_num,int *cuda_real_bnd,int *cuda_where,int *cuda_xadj,\
  int *cuda_adjncy,int *cuda_adjwgt,int *cuda_nparts)
  {
    int ii,iii,pi,id,ed,nnbrs,me,other,k,j,i,me_part,other_part;
    int nparts=cuda_nparts[0];
  
    ii=blockIdx.x*blockDim.x+threadIdx.x,iii= blockIdx.x*blockDim.x+threadIdx.x+1;
  
    for(pi=cuda_ppnvtxs[ii];pi<cuda_ppnvtxs[iii];pi++)
    {
       me=0;
       other=0;
      me_part=cuda_where[pi];
       for(i=cuda_xadj[pi];i<cuda_xadj[pi+1];i++)
       {
             if(cuda_where[cuda_adjncy[i]]==me_part)
             me+=cuda_adjwgt[i];
             else
             other+=cuda_adjwgt[i];
       }

       if(other>me)
       {cuda_real_bnd[atomicAdd(&cuda_real_bnd_num[0],1)]=pi;

      }
    }
  }





__global__ void initpwgts(int *cuda_pwgts)
{
  int  i, ii,iii;
  ii=blockIdx.x*blockDim.x+threadIdx.x;
  cuda_pwgts[ii]=0;

}

__global__ void inittpwgts(float *cuda_ctrl_tpwgts,float *temp)
{
  int  i, ii,iii;
  ii=blockIdx.x*blockDim.x+threadIdx.x;
  cuda_ctrl_tpwgts[ii]=temp[0];

}


__global__ void Sumpwgts(int *cuda_pwgts,int *cuda_where,int *cuda_vwgt)
{
  int  i, ii,iii;
  ii=blockIdx.x*blockDim.x+threadIdx.x;
  atomicAdd(&cuda_pwgts[cuda_where[ii]],cuda_vwgt[ii]);
}

void Mallocinit_refineinfo(ctrl_t *ctrl,graph_t *graph)
{
  int nvtxs=graph->nvtxs;
  int nedges=graph->nedges;
  int nparts=ctrl->nparts;

  t_ctog_start=clock(); 
cudaDeviceSynchronize();
  cudaMalloc((void**)&graph->cuda_where, nvtxs*sizeof(int));
  cudaMemcpy(graph->cuda_where,graph->where,nvtxs*sizeof(int),cudaMemcpyHostToDevice);
cudaDeviceSynchronize();
  t_ctog_end=clock(); 
 duration_ctog += (double)(t_ctog_end - t_ctog_start) / CLOCKS_PER_SEC/test_time; 
  cudaMalloc((void**)&graph->cuda_bnd, nvtxs*sizeof(int));

  int num=0;
  cudaMalloc((void**)&graph->cuda_bndnum, sizeof(int));
  cudaMemcpy(graph->cuda_bndnum,&num,sizeof(int),cudaMemcpyHostToDevice);

 
  cudaMalloc((void**)&graph->cuda_nparts, sizeof(int));
  cudaMemcpy(graph->cuda_nparts,&nparts,sizeof(int),cudaMemcpyHostToDevice);

  cudaMalloc((void**)&graph->cuda_pwgts, ctrl->nparts*sizeof(int));
  initpwgts<<<nparts,1>>>(graph->cuda_pwgts);

  Sumpwgts<<<nvtxs,1>>>(graph->cuda_pwgts,graph->cuda_where,graph->cuda_vwgt);
 
  cudaMalloc((void**)&graph->cuda_tvwgt, sizeof(int));
  cudaMemcpy(graph->cuda_tvwgt,graph->tvwgt,sizeof(int),cudaMemcpyHostToDevice);

  cudaMalloc((void**)&graph->cuda_ctrl_tpwgts, nparts*sizeof(float));

  cudaMalloc((void**)&graph->cuda_maxwgt, nparts*sizeof(int));
  cudaMalloc((void**)&graph->cuda_minwgt, nparts*sizeof(int));
  
  float *temp;
  cudaMalloc((void**)&temp, sizeof(float));
  cudaMemcpy(temp,ctrl->tpwgts,sizeof(int),cudaMemcpyHostToDevice);

  cudaMalloc((void**)&graph->cuda_ctrl_tpwgts, nparts*sizeof(float));

  inittpwgts<<<nparts,1>>>(graph->cuda_ctrl_tpwgts,temp);

  cudaFree(temp);

}

void Mallocrefineinfo (ctrl_t *ctrl,graph_t *graph)
{
  int nvtxs=graph->nvtxs;
  int nedges=graph->nedges;
  int nparts=ctrl->nparts;
  cudaMalloc((void**)&graph->cuda_bnd, nvtxs*sizeof(int));

  
  int num=0;
  cudaMalloc((void**)&graph->cuda_bndnum, sizeof(int));
  cudaMemcpy(graph->cuda_bndnum,&num,sizeof(int),cudaMemcpyHostToDevice);
  
   
  cudaMalloc((void**)&graph->cuda_nparts, sizeof(int));
  cudaMemcpy(graph->cuda_nparts,&nparts,sizeof(int),cudaMemcpyHostToDevice);
  
  cudaMalloc((void**)&graph->cuda_pwgts, ctrl->nparts*sizeof(int));
  initpwgts<<<nparts,1>>>(graph->cuda_pwgts);
  
  Sumpwgts<<<nvtxs,1>>>(graph->cuda_pwgts,graph->cuda_where,graph->cuda_vwgt);

  cudaMalloc((void**)&graph->cuda_tvwgt, sizeof(int));
  cudaMemcpy(graph->cuda_tvwgt,graph->tvwgt,sizeof(int),cudaMemcpyHostToDevice);
  
  cudaMalloc((void**)&graph->cuda_ctrl_tpwgts, nparts*sizeof(float));
  
  cudaMalloc((void**)&graph->cuda_maxwgt, nparts*sizeof(int));
  cudaMalloc((void**)&graph->cuda_minwgt, nparts*sizeof(int));
    
  float *temp;
  cudaMalloc((void**)&temp, sizeof(float));
  cudaMemcpy(temp,ctrl->tpwgts,sizeof(int),cudaMemcpyHostToDevice);
  
  cudaMalloc((void**)&graph->cuda_ctrl_tpwgts, nparts*sizeof(float));
  
  inittpwgts<<<nparts,1>>>(graph->cuda_ctrl_tpwgts,temp);
  cudaFree(temp);
  
}



__global__ void Find_real_bnd_info(int *cuda_real_bnd_num,int *cuda_real_bnd,int *cuda_where,int *cuda_xadj,\
  int *cuda_adjncy,int *cuda_adjwgt,int *cuda_nparts)
  {
    int ii,iii,pi,id,ed,nnbrs,me,other,k,j,i,me_part,other_part;
    int nparts=cuda_nparts[0];
    pi=blockIdx.x*blockDim.x+threadIdx.x;
  
      me=0;
      other=0;
      me_part=cuda_where[pi];
       for(i=cuda_xadj[pi];i<cuda_xadj[pi+1];i++)
       {
             if(cuda_where[cuda_adjncy[i]]==me_part)
             me+=cuda_adjwgt[i];
             else
             other+=cuda_adjwgt[i];
       }

       if(other>me)
       {
         cuda_real_bnd[atomicAdd(&cuda_real_bnd_num[0],1)]=pi;
      }
    
  }

  __global__ void Find_real_bnd_info1(int *cuda_real_bnd_num,int *cuda_real_bnd,int *cuda_where,int *cuda_xadj,\
    int *cuda_adjncy,int *cuda_adjwgt,int *cuda_nparts)
    {
      int ii,iii,pi,id,ed,nnbrs,me,other,k,j,i,me_part,other_part,pii;
      int nparts=cuda_nparts[0];
      pii=blockIdx.x*blockDim.x+threadIdx.x;


        pi=cuda_real_bnd[pii];
        me=0;
        other=0;
        me_part=cuda_where[pi];
         for(i=cuda_xadj[pi];i<cuda_xadj[pi+1];i++)
         {
               if(cuda_where[cuda_adjncy[i]]==me_part)
               me+=cuda_adjwgt[i];
               else
               other+=cuda_adjwgt[i];
         }
  
         if(other>me)
         {
           cuda_real_bnd[atomicAdd(&cuda_real_bnd_num[0],1)]=pi;
        }
  
    }

//ywt
  __global__ void  find_kayparams(int *cuda_real_bnd_num,int *bnd_info,int *cuda_real_bnd,int *cuda_where,int *cuda_xadj,\
    int *cuda_adjncy,int *cuda_adjwgt,int *cuda_nparts,int *cuda_bn,int *cuda_bt,int *cuda_g)
    {
      int ii,iii,pi,id,ed,nnbrs,me,other,k,j,i,me_wgt,other_wgt;
      int start,end;
      int nparts=cuda_nparts[0];
    
      ii=blockIdx.x*blockDim.x+threadIdx.x,iii= blockIdx.x*blockDim.x+threadIdx.x+1;
  
      pi=cuda_real_bnd[ii];

      start=(cuda_nparts[0])*ii;
      end=(cuda_nparts[0])*(ii+1);
  
      for(i=start;i<end;i++)
      bnd_info[i]=0;
  
      for(i=cuda_xadj[pi];i<cuda_xadj[pi+1];i++)
      {
         bnd_info[start+cuda_where[cuda_adjncy[i]]]+=cuda_adjwgt[i];
      }
  
      me_wgt=other_wgt=bnd_info[start+cuda_where[pi]];
      other=cuda_where[pi];
      for(i=start;i<end;i++)
      {
        if(bnd_info[i]>other_wgt)
        {
          other_wgt=bnd_info[i];
          other=i-start;
        }
      }
      cuda_g[ii]=other_wgt-me_wgt;
      cuda_bt[ii]=other;
      cuda_bn[ii]=pi;
  
  }

  __global__ void initbndnum(int *n)
  {
    n[0]=0;
  }


  int refine_pass=1 ;

  __global__ void findcsr(int *bt,int *n,int *nparts,int *bnd_num)
  {
    int ii,iii,pi,id,ed,nnbrs,me,other,k,j,i,me_wgt,other_wgt;
    int start,end;
    ii=blockIdx.x*blockDim.x+threadIdx.x,iii= blockIdx.x*blockDim.x+threadIdx.x+1;

    n[2*ii]=-1;
    n[ 2*ii+1]=-1;
    int temp;

    for(int i=0;i<bnd_num[0];i++)
    {
           if(ii==bt[i])
           {
             n[2*ii]=i;
          break; 
          }

    }
    if(n[2*ii]!=-1)
    {
      for(int i=n[2*ii];i<bnd_num[0];i++)
    {
      if(bt[i]!=ii)
      {
        n[2*ii+1]=i-1;
        break; 
     }
    }
  }
  n[2*bt[bnd_num[0]-1]+1]=bnd_num[0]-1;


  }

  __global__ void check1(int *bnd_num,int *bt)
  {
    for(int i=0;i<bnd_num[0];i++)
    {
      printf("i=%d %d\n",i,bt[i]);
    }
  }


  __global__ void check11(int *bt)
  {
     for(int i=0;i<16;i++)
     printf("csr=%d ",bt[i]);

     printf("\n");
  }

  __global__ void initcucsr(int *cu_csr,int *bndnum)
  {
    cu_csr[0]=0;
    cu_csr[1]=bndnum[0];
  }

__device__ void quicksort_initial(int *arr,int low ,int high,int *arr2)
{
   if(low<high)
   {
      int i=low;
      int j=high;
      int k=arr[low];
      int k2=arr2[low];
      while(i<j)
      {
          while(i<j&&arr[j]>=k)
          {
            j--;
          }
          if(i<j)
          {
            arr[i]=arr[j];
            arr2[i]=arr2[j];
            i++;
          }
          while(i<j&&arr[i]<k)
          {
            i++;
          }
          if(i<j)
          {
            arr[j]=arr[i];
            arr2[j]=arr2[i];
            j--;
          }
      }
      arr[i]=k;
      arr2[i]=k2;
      quicksort_initial(arr,low,i-1,arr2);
      quicksort_initial(arr,i+1,high,arr2);
   }
}

__global__ void sort(int *cu_bt, int *cu_bn,int *bnd_num,int *cu_csr)
  {
     quicksort_initial(cu_bt,0,bnd_num[0]-1,cu_bn);

  }


void Findgraphbndinfo(ctrl_t *ctrl,graph_t *graph)
{


 int nvtxs=graph->nvtxs;
 int nparts=ctrl->nparts;

  cudaDeviceSynchronize();
 initbndnum<<<1,1>>>(graph->cuda_bndnum);
 cudaDeviceSynchronize();

 t_rrk1_start=clock();

 cudaDeviceSynchronize();
 Find_real_bnd_info<<<nvtxs,1>>>(graph->cuda_bndnum,graph->cuda_bnd,graph->cuda_where,graph->cuda_xadj,\
   graph->cuda_adjncy,graph->cuda_adjwgt,graph->cuda_nparts);
   int bnd_num;
   cudaDeviceSynchronize();

   t_rrk1_end=clock();
   duration_rrk1 += (double)(t_rrk1_end - t_rrk1_start) / CLOCKS_PER_SEC/test_time; 

   cudaMemcpy(   &bnd_num,graph->cuda_bndnum, sizeof(int), cudaMemcpyDeviceToHost);

   t_rrk2_start=clock();
   cudaDeviceSynchronize();

   if(bnd_num>0)
   {
   find_kayparams<<<bnd_num,1>>>(graph->cuda_bndnum,cu_info,graph->cuda_bnd,graph->cuda_where,\
     graph->cuda_xadj,graph->cuda_adjncy,graph->cuda_adjwgt,graph->cuda_nparts,\
     cu_bn,cu_bt,cu_g);
   
     cudaDeviceSynchronize();
     t_rrk2_end=clock();
     duration_rrk2 += (double)(t_rrk2_end - t_rrk2_start) / CLOCKS_PER_SEC/test_time; 
    initcucsr<<<1,1>>>(cu_csr,graph->cuda_bndnum);

     bb_segsort(cu_bt, cu_bn,bnd_num,cu_csr,  1);
  //   sort<<<1,1>>>(cu_bt, cu_bn,graph->cuda_bndnum,cu_csr) ;
     findcsr<<<nparts,1>>>(cu_bt,cu_que,graph->cuda_nparts,graph->cuda_bndnum);
  
}
     graph->cpu_bndnum=(int *)malloc(sizeof(int));
     graph->cpu_bndnum[0]=bnd_num;



}


__global__ void Sum_maxmin_pwgts(int *cuda_maxwgt,int *cuda_minwgt,float *cuda_ctrl_tpwgts,int *cuda_tvwgt)
{
  int  i, ii,iii;
  ii=blockIdx.x*blockDim.x+threadIdx.x;

  float ubfactor=1.03;

  cuda_maxwgt[ii]  = int(cuda_ctrl_tpwgts[ii]*cuda_tvwgt[0]*ubfactor);
  cuda_minwgt[ii]  = int(cuda_ctrl_tpwgts[ii]*cuda_tvwgt[0]/ ubfactor);

}


__global__ void Sum_maxmin_pwgts1(int *cuda_maxwgt,int *cuda_minwgt,float *cuda_ctrl_tpwgts,int *cuda_tvwgt)
{
  int  i, ii,iii;
  ii=blockIdx.x*blockDim.x+threadIdx.x;

  float ubfactor=1.03;

  cuda_maxwgt[ii]  = int(cuda_ctrl_tpwgts[ii]*cuda_tvwgt[0]*ubfactor);
  cuda_minwgt[ii]  = int(cuda_ctrl_tpwgts[ii]*cuda_tvwgt[0]/ ubfactor);

}

__global__ void Exchangebndnode(int *pwgts,int *bndnum,int *bnd,int *bndto,int *vwgt,\
  int *maxvwgt,int *minvwgt,int *where)
{
  int me,to,vvwgt,memax,memin,tomax,tomin;
  
  int nmoves=0;

  for(int i=0;i<bndnum[0];i++)
  {
    
    vvwgt=vwgt[bnd[i]];
    me=where[bnd[i]];
    to=bndto[i];
  
    memax=maxvwgt[me];
    memin=minvwgt[me];
    tomax=maxvwgt[to];
    tomin=minvwgt[to];

  
    if(((pwgts[to]+vvwgt>=tomin)&&(pwgts[to]+vvwgt<=tomax))\
    &&((pwgts[me]-vvwgt>=memin)&&(pwgts[me]-vvwgt<=memax)))
  
    {
      
      pwgts[me]-=vvwgt;
      pwgts[to]+=vvwgt;
      where[bnd[i]]=to;
      nmoves++;
  
   }

    else    
   continue;
   
  }


}





__global__ void Exnode1(int *que,int *pwgts,int *bndnum,int *bnd,int *bndto,int *vwgt,\
  int *maxvwgt,int *minvwgt,int *where)
{
  int me,to,vvwgt,memax,memin,tomax,tomin;
  
  int nmoves=0;
  int  i, ii,iii;
  ii=blockIdx.x*blockDim.x+threadIdx.x;

  if(que[2*ii]!=-1)
  {
  for(i=que[2*ii];i<=que[2*ii+1];i++)
  {
 
    vvwgt=vwgt[bnd[i]];
    me=where[bnd[i]];
    to=bndto[i];

//    printf("%d %d\n",me ,to);
  
    memax=maxvwgt[me];
    memin=minvwgt[me];
    tomax=maxvwgt[to];
    tomin=minvwgt[to];

  if(me<=to)
  {
    if(((pwgts[to]+vvwgt>=tomin)&&(pwgts[to]+vvwgt<=tomax))\
    &&((pwgts[me]-vvwgt>=memin)&&(pwgts[me]-vvwgt<=memax)))
    {
 
      atomicAdd(&pwgts[to],vvwgt);
      atomicSub(&pwgts[me],vvwgt);

      where[bnd[i]]=to;
      nmoves++;
  
   }
  
  }
  
 }
}

}

__global__ void Exnode2(int *que,int *pwgts,int *bndnum,int *bnd,int *bndto,int *vwgt,\
  int *maxvwgt,int *minvwgt,int *where)
{
  int me,to,vvwgt,memax,memin,tomax,tomin;
  
  int nmoves=0;
  int  i, ii,iii;
  ii=blockIdx.x*blockDim.x+threadIdx.x;

  if(que[2*ii]!=-1)
  {
  for(i=que[2*ii];i<=que[2*ii+1];i++)
  {
    
    vvwgt=vwgt[bnd[i]];
    me=where[bnd[i]];
    to=bndto[i];
  
    memax=maxvwgt[me];
    memin=minvwgt[me];
    tomax=maxvwgt[to];
    tomin=minvwgt[to];

  if(me>to)
  {
    if(((pwgts[to]+vvwgt>=tomin)&&(pwgts[to]+vvwgt<=tomax))\
    &&((pwgts[me]-vvwgt>=memin)&&(pwgts[me]-vvwgt<=memax)))
  
    {
      
      atomicAdd(&pwgts[to],vvwgt);
      atomicSub(&pwgts[me],vvwgt);
      where[bnd[i]]=to;
      nmoves++;
  
   }
  
  }
  }
  }
 
}



void Kwayrefine(ctrl_t *ctrl,graph_t *graph)
{
 
  int nvtxs=graph->nvtxs; 
  int nparts=ctrl->nparts;
  cudaDeviceSynchronize();
  Sum_maxmin_pwgts<<<nparts,1>>>(graph->cuda_maxwgt,graph->cuda_minwgt,graph->cuda_ctrl_tpwgts,\
    graph->cuda_tvwgt);

 
    cudaDeviceSynchronize();

    for(int i=0;i<refine_pass;i++)
    {
       Findgraphbndinfo(ctrl,graph);

       if(graph->cpu_bndnum[0]>0)
      {
        Exnode1<<<ctrl->nparts,1>>>(cu_que,graph->cuda_pwgts,graph->cuda_bndnum,cu_bn,cu_bt,graph->cuda_vwgt,\
          graph->cuda_maxwgt,graph->cuda_minwgt,graph->cuda_where);
          
          Exnode2<<<ctrl->nparts,1>>>(cu_que,graph->cuda_pwgts,graph->cuda_bndnum,cu_bn,cu_bt,graph->cuda_vwgt,\
            graph->cuda_maxwgt,graph->cuda_minwgt,graph->cuda_where);
        
      }
      else
      break;


    }



}



void Kwayrefine1(ctrl_t *ctrl,graph_t *graph)
{
 
  int nvtxs=graph->nvtxs; 
  int nparts=ctrl->nparts;
  cudaDeviceSynchronize();
  Sum_maxmin_pwgts1<<<nparts,1>>>(graph->cuda_maxwgt,graph->cuda_minwgt,graph->cuda_ctrl_tpwgts,\
    graph->cuda_tvwgt);
    cudaDeviceSynchronize();

    for(int i=0;i<1;i++)
    {
       Findgraphbndinfo(ctrl,graph);
 //      printf("bnd=%d\n",graph->cpu_bndnum[0]);

       if(graph->cpu_bndnum[0]>0)
      {
        Exnode1<<<ctrl->nparts,1>>>(cu_que,graph->cuda_pwgts,graph->cuda_bndnum,cu_bn,cu_bt,graph->cuda_vwgt,\
          graph->cuda_maxwgt,graph->cuda_minwgt,graph->cuda_where);
          Findgraphbndinfo(ctrl,graph);
          if(graph->cpu_bndnum[0]>0)
          {
          Exnode2<<<ctrl->nparts,1>>>(cu_que,graph->cuda_pwgts,graph->cuda_bndnum,cu_bn,cu_bt,graph->cuda_vwgt,\
            graph->cuda_maxwgt,graph->cuda_minwgt,graph->cuda_where);
          }
      }
      else
      break;


    }



}


void Freecudarefine(graph_t *graph)//ywt
{
  cudaFree(graph->cuda_real_cnedges);
  cudaFree(graph->cuda_xadj);
  cudaFree(graph->cuda_cmap);
  cudaFree(graph->cuda_nvtxs);
  cudaFree(graph->cuda_adjncy);
  cudaFree(graph->cuda_adjwgt);
  cudaFree(graph->cuda_vwgt);

  cudaFree(graph->cuda_maxwgt);
  cudaFree(graph->cuda_minwgt);
  cudaFree(graph->cuda_where);
  cudaFree(graph->cuda_pwgts);
  cudaFree(graph->cuda_bnd);
  cudaFree(graph->cuda_bndto);
  cudaFree(graph->cuda_bndnum);
  //cudaFree(graph->cuda_gain);
  cudaFree(graph->bnd_info);
  cudaFree(graph->cuda_real_bnd_num);
  cudaFree(graph->cuda_real_bnd);
  cudaFree(graph->cuda_nparts);
  cudaFree(graph->cuda_bn);
  cudaFree(graph->cuda_bt);
  cudaFree(graph->cuda_g);
  cudaFree(graph->cuda_tvwgt);
  cudaFree(graph->cuda_ctrl_tpwgts);

}

// __global__ void check(int *a)
// {
//   printf("check=%d ",a[0]);
// }

// void test426(graph_t *graph)
// {
//   printf("%d ",graph->nvtxs);
//   //cudaFree(graph->cuda_maxwgt);
//   cudaDeviceSynchronize();
//   check<<<1,1>>>(graph->cuda_xadj);
//   cudaDeviceSynchronize();

// }





void InitialPartitioninh(ctrl_t *ctrl,graph_t *graph)
{
  graph->where=(int *)malloc(sizeof(int)*graph->nvtxs);


 int i,j,k;
  int nvtxs=graph->nvtxs;
  int nparts=ctrl->nparts;

  int  *p=(int*)malloc(sizeof(int)*nparts);
  int  *pp=(int*)malloc(sizeof(int)*nparts);
  int  *w=(int*)malloc(sizeof(int)*nparts);
  int  *g=(int*)malloc(sizeof(int)*nparts);
  
  double init_ubfactor=1.03;

  int pwgts=0;
  for(i=0;i<nvtxs;i++)
  pwgts+=graph->vwgt[i];

//  int C=int(ctrl->tpwgts[0]*pwgts/init_ubfactor);

//int C=int(pwgts*init_ubfactor)/nparts;

//int C=int(pwgts/1.001)/nparts;

int C=pwgts/nparts;

printf("C=%d\n",C);

  for(i=0;i<nparts;i++)
  {
    p[i]=0;
    pp[i]=0;
    w[i]=0;
    g[i]=0;
  }


  for(i=0;i<nvtxs;i++)
  {
    graph->where[i]=-1;
  }

graph->where[0]=0; 

  for(i=1;i<nvtxs;i++)
  {
    for(j=graph->xadj[i+1]-1;j>=graph->xadj[i];j--)
     {
       k=graph->adjncy[j];
        if(graph->where[k]!=-1)
        {
            pp[graph->where[k]]+=graph->vwgt[k];
        }
     }

     for(j=0;j<nparts;j++)
    {
      w[j]=C-p[j];
    }
     
    for(j=0;j<nparts;j++)
    {
      if(pp[j]==0&&w[j]<0)
      g[j]=-1;
      else
      g[j]=pp[j]*w[j];
      
    }

    int temp=-1;
    double temp1=-1;
    for(j=0;j<nparts;j++)
    {
        if(g[j]>temp1)
        {
          temp1=g[j];
          temp=j;
        }
    }


    graph->where[i]=temp;
    p[graph->where[i]]+=graph->vwgt[i];

    for(j=0;j<nparts;j++)
    pp[j]=0;

  }

for(i=0;i<nparts;i++)
printf("%d %d\n",i,p[i]);

}

__global__ void checkp(int *n,int *w)
{
  printf("n=%d\n",n[0]);

  
  for(int i=0;i<n[0];i++)
  printf("%d",w[i]);
}

int MlevelKWayPartitioning(ctrl_t *ctrl, graph_t *graph, int *part)
{
  int i, j, objval=0, curobj=0, bestobj=0;
  float curbal=0.0, bestbal=0.0;

  graph_t *cgraph;
  t_coarsen_start=clock(); 

  cgraph = CoarsenGraph(ctrl, graph);

  t_coarsen_end=clock(); 
  duration_coarsen += (double)(t_coarsen_end - t_coarsen_start) / CLOCKS_PER_SEC/test_time; 

 // printf("%d \n",cgraph->nvtxs);
  t_init_start=clock(); 
  printf("init start\n");
  InitKWayPartitioning(ctrl, cgraph);
  printf("init end\n");

  printf("%d",cgraph->nvtxs);
 //InitialPartitioninh(ctrl,cgraph);
 
// int e=ComputeCut(cgraph, cgraph->where);
// printf("edge-cut=%d\n",e);
 

  t_init_end=clock(); 
  duration_init += (double)(t_init_end - t_init_start) / CLOCKS_PER_SEC/test_time; 


  t_uncoarsen_start=clock(); 

  // for(i=0;i<cgraph->nvtxs;i++)
  // printf("%d ",cgraph->where[i]);

  // printf("\n\n\n\n\n\n");

  
  Mallocinit_refineinfo(ctrl,cgraph);

//Kwayrefine1(ctrl,cgraph);

//cudaMemcpy(   cgraph->where,cgraph->cuda_where, cgraph->nvtxs*sizeof(int), cudaMemcpyDeviceToHost);

//e=ComputeCut(cgraph, cgraph->where);
// printf("edge-cut=%d\n",e);
/*
  cudaDeviceSynchronize();
  checkp<<<1,1>>>(cgraph->cuda_nvtxs,cgraph->cuda_where);  
  cudaDeviceSynchronize();
*/
  Kwayrefine(ctrl,cgraph);

  for(i=0;;i++)
  {
  
  if(cgraph!=graph)
  {
//Freecudarefine(cgraph);
    cgraph=cgraph->finer;
    cudaMalloc((void**)&cgraph->cuda_where, cgraph->nvtxs*sizeof(int));
    ProjectKWayPartition(ctrl,cgraph); 
    Mallocrefineinfo(ctrl,cgraph);
    Kwayrefine(ctrl,cgraph);
    Freecudarefine(cgraph->coarser);

    //test426(cgraph->coarser);

 // printf("\n");
  }
  else 
  {
  
  break; 
  }
  }
  t_uncoarsen_end=clock(); 
  duration_uncoarsen += (double)(t_uncoarsen_end - t_uncoarsen_start) / CLOCKS_PER_SEC/test_time; 
  
  return 0;
}





// int PartGraph(int *nvtxs,  int *xadj, int *adjncy, int *vwgt,int *adjwgt, \
// int *nparts, float *tpwgts, float *ubvec, int *objval, int *part)
// {
          
//     graph_t *graph;
//     ctrl_t *ctrl;
//     ctrl = SetupCtrl( *nparts, tpwgts, ubvec);
//     graph = SetupGraph(ctrl, *nvtxs, xadj, adjncy, vwgt, adjwgt);
//     SetupKWayBalMultipliers(ctrl, graph);
//     ctrl->CoarsenTo = gk_max((*nvtxs)/(20*gk_log2(*nparts)), 30*(*nparts));



//     ctrl->nIparts   = (ctrl->CoarsenTo == 30*(*nparts) ? 4 : 5);
//     AllocateWorkSpace(ctrl, graph);

//     test_time=1;


//     Mallocoriginal_coarseninfo(ctrl,graph);
    
//      cudaMalloc((void**)&cu_info, graph->nvtxs*ctrl->nparts* sizeof(int));
//      cudaMalloc((void**)&cu_bn, graph->nvtxs* sizeof(int));
//      cudaMalloc((void**)&cu_bt, graph->nvtxs* sizeof(int));
//      cudaMalloc((void**)&cu_g, graph->nvtxs*sizeof(int));
//      cudaMalloc((void**)&cu_csr, 2*sizeof(int));
//      cudaMalloc((void**)&cu_que, 2*ctrl->nparts*sizeof(int));
     
//     cudaDeviceSynchronize();
//     //time
//     t_part_start=clock();  
//     cudaDeviceSynchronize();

//     for(int i=0;i<test_time;i++)
//     {
//     *objval = MlevelKWayPartitioning(ctrl, graph, part);
//     }
//     cudaDeviceSynchronize();
//     t_part_end=clock();
//     duration_part = (double)(t_part_end - t_part_start) / CLOCKS_PER_SEC/test_time; 
    
//     cudaMemcpy(   part,graph->cuda_where, graph->nvtxs*sizeof(int), cudaMemcpyDeviceToHost);
//     cudaDeviceSynchronize();
    
//     Freecudacoarsen(graph);
//     Freecudarefine(graph);
    
//   cudaFree(cu_info);
//   cudaFree(cu_bn);
//   cudaFree(cu_bt);
//   cudaFree(cu_g);
//   cudaFree(cu_csr);
//   cudaFree(cu_que);

//   //test426(graph);

 
//     return 1;
// }




// void errexit(char *f_str,...)
// {
//   va_list argp;

//   va_start(argp, f_str);
//   vfprintf(stderr, f_str, argp);
//   va_end(argp);

//   if (strlen(f_str) == 0 || f_str[strlen(f_str)-1] != '\n')
//         fprintf(stderr,"\n");
//   fflush(stderr);

//   if (1)
//     exit(-2);

// }

// FILE *gk_fopen(char *fname, char *mode, const char *msg)
// {
//   FILE *fp;
//   char errmsg[8192];

//   fp = fopen(fname, mode);
//   if (fp != NULL)
//     return fp;

//   sprintf(errmsg,"file: %s, mode: %s, [%s]", fname, mode, msg);
//   perror(errmsg);
//   errexit("Failed on gk_fopen()\n");

//   return NULL;
// }

// int gk_getline(char **lineptr, size_t *n, FILE *stream)
// {
//   return getline(lineptr, n, stream);
// }
// graph_t *ReadGraph(char *filename)
// {
//   int i, j, k, l, fmt, nfields, readew, readvw, readvs, edge, ewgt;
//   int *xadj, *adjncy, *vwgt, *adjwgt, *vsize;
//   char *line=NULL, fmtstr[256], *curstr, *newstr;
//   size_t lnlen=0;
//   FILE *fpin;
//   graph_t *graph;

//   graph = CreateGraph();
    
//   fpin = gk_fopen(filename, "r", "ReadGRaph: Graph");


//      do {
//     if (gk_getline(&line, &lnlen, fpin) == -1) 
//       errexit("Premature end of input file: file: %s\n", filename);

//   } while (line[0] == '%');
     

//   fmt= 0;
  
//   nfields = sscanf(line, "%d %d %d", &(graph->nvtxs), &(graph->nedges), &fmt);
  

//   if (nfields < 2) 
//     errexit("The input file does not specify the number of vertices and edges.\n");

//   if (graph->nvtxs <= 0 || graph->nedges <= 0) 
//    errexit("The supplied nvtxs:%d and nedges:%d must be positive.\n",  graph->nvtxs, graph->nedges);
        
//   if (fmt > 111) 
//     errexit("Cannot read this type of file format [fmt=%d]!\n", fmt);

//   sprintf(fmtstr, "%03d", fmt%1000);
//   readvs = (fmtstr[0] == '1');
//   readvw = (fmtstr[1] == '1');
//   readew = (fmtstr[2] == '1');
    



//   graph->nedges *=2;

//   xadj   = graph->xadj   = (int*)malloc(sizeof(int)*(graph->nvtxs+1));
//   for(i=0;i<graph->nvtxs+1;i++)
//   {
//     xadj[i]=graph->xadj[i]=0;
//   }
 
//   adjncy = graph->adjncy=(int*)malloc(sizeof(int)*(graph->nedges));
  

//   vwgt   = graph->vwgt   =  (int*)malloc(sizeof(int)*(graph->nvtxs));
//     for(i=0;i<graph->nvtxs;i++)
//   {
//     vwgt[i]=graph->vwgt[i]=1;
//   }

//   adjwgt = graph->adjwgt = (int*)malloc(sizeof(int)*(graph->nedges));
//       for(i=0;i<graph->nedges;i++)
//   {
//     adjwgt[i]=graph->adjwgt[i]=1;
//   }



//   for (xadj[0]=0, k=0, i=0; i<graph->nvtxs; i++) {
 
//     do {
      
//       if (getline(&line, &lnlen, fpin) == -1) 
//   errexit("Premature end of input file while reading vertex %d.\n", i+1); 
//     } while (line[0] == '%');

//     curstr = line;
//     newstr = NULL;

  


//     if (readvw) {

//         vwgt[i] = strtol(curstr, &newstr, 10);
//         if (newstr == curstr)
//           errexit("The line for vertex %d does not have enough weights "
//                   "for the %d constraints.\n", i+1, 1);
//         if (vwgt[i] < 0)
//           errexit("The weight vertex %d and constraint %d must be >= 0\n", i+1, 0);
//         curstr = newstr;
      
//     }
 
//     while (1) {
//       edge = strtol(curstr, &newstr, 10);
//       if (newstr == curstr)
//         break; 
//       curstr = newstr;

//       if (edge < 1 || edge > graph->nvtxs)
//         errexit("Edge %d for vertex %d is out of bounds\n", edge, i+1);

//       ewgt = 1;
//       if (readew) {
//         ewgt = strtol(curstr, &newstr, 10);
//         if (newstr == curstr)
//           errexit("Premature end of line for vertex %d\n", i+1);
//         if (ewgt <= 0)
//           errexit("The weight (%d) for edge (%d, %d) must be positive.\n",    ewgt, i+1, edge);
//         curstr = newstr;
//       }

//       if (k == graph->nedges)
//         errexit("There are more edges in the file than the %d specified.\n", graph->nedges/2);

//       adjncy[k] = edge-1;
//       adjwgt[k] = ewgt;
//       k++;
//     } 

//     xadj[i+1] = k;
//   }
  
//   fclose(fpin);

//   if (k != graph->nedges) {
//       printf("------------------------------------------------------------------------------\n");
//     printf("***  I detected an error in your input file  ***\n\n");
//     printf("In the first line of the file, you specified that the graph contained\n"
//            "%d edges. However, I only found %d edges in the file.\n", 
//            graph->nedges/2, k/2);
//     if (2*k == graph->nedges) {
//       printf("\n *> I detected that you specified twice the number of edges that you have in\n");
//       printf("    the file. Remember that the number of edges specified in the first line\n");
//       printf("    counts each edge between vertices v and u only once.\n\n");
//     }
//     printf("Please specify the correct number of edges in the first line of the file.\n");
//     printf("------------------------------------------------------------------------------\n");
//     exit(0);
  
//   }

//   gk_free((void **)&line, LTERM);

//   return graph;
// }


// void WritePartition(char *fname, int *part, int n, int nparts)
// {
//   FILE *fpout;
//   int i;
//   char filename[1280000];

//   sprintf(filename, "%s.part.%d", fname, nparts);

//   fpout = gk_fopen(filename, "w", __func__);

//   for (i=0; i<n; i++)
//     fprintf(fpout,"%d\n", part[i]);

//   fclose(fpout);
// }



int PartGraph(int *nvtxs,  int *xadj, int *adjncy, int *vwgt,int *adjwgt, \
int *nparts, float *tpwgts, float *ubvec, int *objval, int *part)
{
          
    graph_t *graph;
    ctrl_t *ctrl;
    ctrl = SetupCtrl( *nparts, tpwgts, ubvec);
    graph = SetupGraph(ctrl, *nvtxs, xadj, adjncy, vwgt, adjwgt);
    SetupKWayBalMultipliers(ctrl, graph);
    ctrl->CoarsenTo = gk_max((*nvtxs)/(20*gk_log2(*nparts)), 30*(*nparts));



    ctrl->nIparts   = (ctrl->CoarsenTo == 30*(*nparts) ? 4 : 5);
    AllocateWorkSpace(ctrl, graph);

    test_time=1;


    Mallocoriginal_coarseninfo(ctrl,graph);
    
    cudaMalloc((void**)&cu_info, graph->nvtxs * ctrl->nparts * sizeof(int));
     


    double data_byte =NULL;
    double data_MB = NULL;
    printf("graph->nvtxs = %d,ctrl->nparts = %d\n",graph->nvtxs,ctrl->nparts);
    
    data_byte = (double)graph->nvtxs * (double)ctrl->nparts * 4;
    
    data_MB = data_byte /1024/1024;

    printf("data_byte = %lf byte, data_MB = %lf MB\n",data_byte,data_MB);






     cudaMalloc((void**)&cu_bn, graph->nvtxs* sizeof(int));
     cudaMalloc((void**)&cu_bt, graph->nvtxs* sizeof(int));
     cudaMalloc((void**)&cu_g, graph->nvtxs*sizeof(int));
     cudaMalloc((void**)&cu_csr, 2*sizeof(int));
     cudaMalloc((void**)&cu_que, 2*ctrl->nparts*sizeof(int));
     
    cudaDeviceSynchronize();
    //time
    t_part_start=clock();  
    cudaDeviceSynchronize();

    for(int i=0;i<test_time;i++)
    {
    *objval = MlevelKWayPartitioning(ctrl, graph, part);
    }
    cudaDeviceSynchronize();
    t_part_end=clock();
    duration_part = (double)(t_part_end - t_part_start) / CLOCKS_PER_SEC/test_time; 
    
    cudaMemcpy(   part,graph->cuda_where, graph->nvtxs*sizeof(int), cudaMemcpyDeviceToHost);
    cudaDeviceSynchronize();
    
    Freecudacoarsen(graph);
    Freecudarefine(graph);
    
  cudaFree(cu_info);
  cudaFree(cu_bn);
  cudaFree(cu_bt);
  cudaFree(cu_g);
  cudaFree(cu_csr);
  cudaFree(cu_que);

  //test426(graph);

 
    return 1;
}




void errexit(char *f_str,...)
{
  va_list argp;

  va_start(argp, f_str);
  vfprintf(stderr, f_str, argp);
  va_end(argp);

  if (strlen(f_str) == 0 || f_str[strlen(f_str)-1] != '\n')
        fprintf(stderr,"\n");
  fflush(stderr);

  if (1)
    exit(-2);

}

FILE *gk_fopen(char *fname, char *mode, const char *msg)
{
  FILE *fp;
  char errmsg[8192];

  fp = fopen(fname, mode);
  if (fp != NULL)
    return fp;

  sprintf(errmsg,"file: %s, mode: %s, [%s]", fname, mode, msg);
  perror(errmsg);
  errexit("Failed on gk_fopen()\n");

  return NULL;
}

int gk_getline(char **lineptr, size_t *n, FILE *stream)
{
  return getline(lineptr, n, stream);
}
graph_t *ReadGraph(char *filename)
{
  int i, j, k, l, fmt, nfields, readew, readvw, readvs, edge, ewgt;
  int *xadj, *adjncy, *vwgt, *adjwgt, *vsize;
  char *line=NULL, fmtstr[256], *curstr, *newstr;
  size_t lnlen=0;
  FILE *fpin;
  graph_t *graph;

  graph = CreateGraph();
    
  fpin = gk_fopen(filename, "r", "ReadGRaph: Graph");


     do {
    if (gk_getline(&line, &lnlen, fpin) == -1) 
      errexit("Premature end of input file: file: %s\n", filename);

  } while (line[0] == '%');
     

  fmt= 0;
  
  nfields = sscanf(line, "%d %d %d", &(graph->nvtxs), &(graph->nedges), &fmt);
  

  if (nfields < 2) 
    errexit("The input file does not specify the number of vertices and edges.\n");

  if (graph->nvtxs <= 0 || graph->nedges <= 0) 
   errexit("The supplied nvtxs:%d and nedges:%d must be positive.\n",  graph->nvtxs, graph->nedges);
        
  if (fmt > 111) 
    errexit("Cannot read this type of file format [fmt=%d]!\n", fmt);

  sprintf(fmtstr, "%03d", fmt%1000);
  readvs = (fmtstr[0] == '1');
  readvw = (fmtstr[1] == '1');
  readew = (fmtstr[2] == '1');
    



  graph->nedges *=2;

  xadj   = graph->xadj   = (int*)malloc(sizeof(int)*(graph->nvtxs+1));
  for(i=0;i<graph->nvtxs+1;i++)
  {
    xadj[i]=graph->xadj[i]=0;
  }
 
  adjncy = graph->adjncy=(int*)malloc(sizeof(int)*(graph->nedges));
  

  vwgt   = graph->vwgt   =  (int*)malloc(sizeof(int)*(graph->nvtxs));
    for(i=0;i<graph->nvtxs;i++)
  {
    vwgt[i]=graph->vwgt[i]=1;
  }

  adjwgt = graph->adjwgt = (int*)malloc(sizeof(int)*(graph->nedges));
      for(i=0;i<graph->nedges;i++)
  {
    adjwgt[i]=graph->adjwgt[i]=1;
  }



  for (xadj[0]=0, k=0, i=0; i<graph->nvtxs; i++) {
 
    do {
      
      if (getline(&line, &lnlen, fpin) == -1) 
  errexit("Premature end of input file while reading vertex %d.\n", i+1); 
    } while (line[0] == '%');

    curstr = line;
    newstr = NULL;

  


    if (readvw) {

        vwgt[i] = strtol(curstr, &newstr, 10);
        if (newstr == curstr)
          errexit("The line for vertex %d does not have enough weights "
                  "for the %d constraints.\n", i+1, 1);
        if (vwgt[i] < 0)
          errexit("The weight vertex %d and constraint %d must be >= 0\n", i+1, 0);
        curstr = newstr;
      
    }
 
    while (1) {
      edge = strtol(curstr, &newstr, 10);
      if (newstr == curstr)
        break; 
      curstr = newstr;

      if (edge < 1 || edge > graph->nvtxs)
        errexit("Edge %d for vertex %d is out of bounds\n", edge, i+1);

      ewgt = 1;
      if (readew) {
        ewgt = strtol(curstr, &newstr, 10);
        if (newstr == curstr)
          errexit("Premature end of line for vertex %d\n", i+1);
        if (ewgt <= 0)
          errexit("The weight (%d) for edge (%d, %d) must be positive.\n",    ewgt, i+1, edge);
        curstr = newstr;
      }

      if (k == graph->nedges)
        errexit("There are more edges in the file than the %d specified.\n", graph->nedges/2);

      adjncy[k] = edge-1;
      adjwgt[k] = ewgt;
      k++;
    } 

    xadj[i+1] = k;
  }
  
  fclose(fpin);

  if (k != graph->nedges) {
      printf("------------------------------------------------------------------------------\n");
    printf("***  I detected an error in your input file  ***\n\n");
    printf("In the first line of the file, you specified that the graph contained\n"
           "%d edges. However, I only found %d edges in the file.\n", 
           graph->nedges/2, k/2);
    if (2*k == graph->nedges) {
      printf("\n *> I detected that you specified twice the number of edges that you have in\n");
      printf("    the file. Remember that the number of edges specified in the first line\n");
      printf("    counts each edge between vertices v and u only once.\n\n");
    }
    printf("Please specify the correct number of edges in the first line of the file.\n");
    printf("------------------------------------------------------------------------------\n");
    exit(0);
  
  }

  gk_free((void **)&line, LTERM);

  return graph;
}


void WritePartition(char *fname, int *part, int n, int nparts)
{
  FILE *fpout;
  int i;
  char filename[1280000];

  sprintf(filename, "%s.part.%d", fname, nparts);

  fpout = gk_fopen(filename, "w", __func__);

  for (i=0; i<n; i++)
    fprintf(fpout,"%d\n", part[i]);

  fclose(fpout);
}





int *mtx_partition(SSS_IMAT *S, int part)
{
    int i,j;

    idx_t rownum = S->num_rows, weights = 1, objval, parts = part;
    idx_t *xadj = (idx_t *)malloc(sizeof(idx_t)*(rownum+1));
    idx_t *reorder = (idx_t *)malloc(sizeof(idx_t)*rownum);


    float *tpwgts = (float *)malloc(parts *sizeof(float));
    for(int i = 0;i<parts;i++)
    {
        tpwgts[i] = 1 / parts ;
    }
    float ubvec = 1.03;


    xadj[0] = 0;
    for ( i = 0; i < rownum; i++)
    {
        int flag = 0;
        for( j = S->row_ptr[i]; j < S->row_ptr[i+1]; j++)
        {
            if(S->col_idx[j] == i)
            {
                flag = 1;
                break;
            }
        }
        if(flag == 1)
            xadj[i+1] = S->row_ptr[i+1] - S->row_ptr[i] - 1;
        else
            xadj[i+1] = S->row_ptr[i+1] - S->row_ptr[i];
    }




    for (i = 1; i <= rownum; i++) xadj[i] += xadj[i-1];

    idx_t *adjncy= (idx_t *)malloc(sizeof(idx_t)*xadj[rownum]);
    
    for (i = 0; i < rownum; i++)
    {
        int tt = xadj[i];
        int ttt = 0;
        for(j = S->row_ptr[i]; j < S->row_ptr[i+1]; j++)
        {
            if(S->col_idx[j] == i) continue;
            else
            {
                adjncy[tt+ttt] = S->col_idx[j];
                ttt++;
            }
        }
    }

    

    double metis1 = SSS_get_time();
    
    //kway
    // printf("metis_kway :\n");
    //METIS_PartGraphKway(&rownum, &weights, xadj, adjncy, NULL, NULL, NULL, &parts, NULL, NULL, NULL, &objval, reorder);
    
    
    //recu
    printf("metis_rece :\n");
    METIS_PartGraphRecursive(&rownum, &weights, xadj, adjncy, NULL, NULL, NULL, &parts, NULL, NULL, NULL, &objval, reorder);





    


    double metis2 = SSS_get_time();

    metis_time += metis2-metis1;


       


    int *order = (int *)malloc(sizeof(int)*rownum);
    for (i = 0; i < rownum; i++)
        order[i] = reorder[i];



    // int *test = (int *)malloc(sizeof(int ) * part);
    
    // for(int i=0;i<part;i++)
    // {
    //   test[i] = 0;      
    // }

    // for(int i = 0;i<rownum;i++)
    // {
    //   test[order[i]]++;
    // }

    // for(int i=0;i<part;i++)
    // {
    //   printf("test[%d] = %d\n",i,test[i]);
    // }




    // for (i = 0; i < rownum; i++) 
    // printf("%d %d\n",i,order[i]);  

    free(reorder);
    free(adjncy);   
    free(xadj);
    
    return order;
}


int *mtx_partition_device(SSS_IMAT *S, int part)
{
    int i,j;

    int rownum = S->num_rows, weights = 1, objval, parts = part;
    int *xadj = (int *)malloc(sizeof(int)*(rownum+1));
    int *reorder = (int *)malloc(sizeof(int)*rownum);

    float *tpwgts = (float *)malloc(parts *sizeof(float));
    for(int i = 0;i<parts;i++)
    {
        tpwgts[i] = 1 / parts ;
    }
    float ubvec = 1.03;


    xadj[0] = 0;
    for ( i = 0; i < rownum; i++)
    {
        int flag = 0;
        for( j = S->row_ptr[i]; j < S->row_ptr[i+1]; j++)
        {
            if(S->col_idx[j] == i)
            {
                flag = 1;
                break;
            }
        }
        if(flag == 1)
            xadj[i+1] = S->row_ptr[i+1] - S->row_ptr[i] - 1;
        else
            xadj[i+1] = S->row_ptr[i+1] - S->row_ptr[i];
    }




    for (i = 1; i <= rownum; i++) xadj[i] += xadj[i-1];

    int *adjncy= (int *)malloc(sizeof(int)*xadj[rownum]);
    
    
    for (i = 0; i < rownum; i++)
    {
        int tt = xadj[i];
        int ttt = 0;
        for(j = S->row_ptr[i]; j < S->row_ptr[i+1]; j++)
        {
            if(S->col_idx[j] == i) continue;
            else
            {
                adjncy[tt+ttt] = S->col_idx[j];
                ttt++;
            }
        }
    }

    

    double metis1 = SSS_get_time();
    
    //gpu
    printf("GPU metis\n");
    PartGraph(&rownum,  xadj, adjncy, NULL, NULL, &parts,tpwgts,&ubvec, &objval, reorder);

    double metis2 = SSS_get_time();

    metis_time += metis2-metis1;


       


    int *order = (int *)malloc(sizeof(int)*rownum);
    for (i = 0; i < rownum; i++)
        order[i] = reorder[i];


    
    // int *test = (int *)malloc(sizeof(int ) * part);
  
    // for(int i=0;i<part;i++)
    // {
    //   test[i] = 0;      
    // }

    // for(int i = 0;i<rownum;i++)
    // {
    //   test[order[i]]++;
    // }

    // for(int i=0;i<part;i++)
    // {
    //   printf("test[%d] = %d\n",i,test[i]);
    // }





    free(reorder);
    free(adjncy);   
    free(xadj);
    
    return order;
}



SSS_IVEC SSS_amg_agg_GPU(SSS_MAT *A,  SSS_MAT *P, SSS_IMAT *S,SSS_AMG_PARS *pars)
{
    cudaError_t EE;

    SSS_IVEC aggregates;
    int maxnagg;
    
    /* ---------------- metis ---------------- */

    // if (S->num_rows  >= 100000)
    // {
    //     maxnagg = 10000;            
    // }
    // else 
    // {
    //     maxnagg = S->num_rows / 10;            
    // }
    maxnagg = S->num_rows / 3;            


    //划分子图数 maxnagg
    printf("maxnagg = %d \n",maxnagg);

    //节点数 S->num_rows 
    int *metis_res = (int*)malloc( S->num_rows * sizeof(int));

    //CPU_metis
    metis_res = mtx_partition(S,maxnagg);

    //GPU_metis
    //metis_res = mtx_partition_device(S,maxnagg);


    pars->n_agg = maxnagg;
    printf("maxnagg = %d\n",maxnagg);



    aggregates = SSS_ivec_create_device(S->num_rows,aggregates);

    EE = cudaMemcpy(aggregates.d_d,metis_res,S->num_rows * sizeof(int),cudaMemcpyHostToDevice);

    if(EE != cudaSuccess)
    {
        printf("aggregates.d_d fail \n");
    }


    

    return aggregates;
}
